Semiclassical trace formulas for noninteracting identical particles 
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We extend the Gutzwiller trace formula to systems of noninteracting identical particles. The 
standard relation for isolated orbits does not apply since the energy of each particle is separately 
conserved causing the periodic orbits to occur in continuous families. The identical nature of the 
particles also introduces discrete permutational symmetries. We exploit the formalism of Creagh and 
Littlejohn [Phys. Rev. A 44, 836 (1991)], who have studied semiclassical dynamics in the presence 
of continuous symmetries, to derive many-body trace formulas for the full and symmetry-reduced 
densities of states. Numerical studies of the three-particle cardioid billiard are used to explicitly 
illustrate and test the results of the theory. 
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I. INTRODUCTION 

In the semiclassical limit of quantum mechanics, the 
periodic orbits of the corresponding classical system play 
a special role in determining the spectral properties of 
the quantum system. This fundamental fact has been a 
dominant theme in modern semiclassical physics and was 
pioneered by Gutzwiller Q , Balian and Bloch [p| , Struti- 
nsky and Magner || and Berry and Tabor j| . One of the 
central results which emerged from this work is the repre- 
sentation of the density of states in terms of classical peri- 
odic orbits. Such representations are referred to as trace 
formulas. Semiclassical analysis based on the use of trace 
formulas is now common in many areas of physics M, Q . 
Besides providing a natural framework for studying the 
quantum manifestations of classical chaos [§, PL such 
analysis has been used in the study of nuclei pj, jlQ, 11 1, 
atoms [[[2], ^3), metal clusters jl4|, ^5|, molecules [16 , 
chemical systems [jl7|, spins [[l8), Casimir effects and 
tunneling pCfl . Trace formulas have also become a promi- 
nent analytical tool in the study of mesoscopic systems 
plj . New directions continue to be explored p2| . 

Despite the vast utility of trace formulas, their use in 
the few-body or many-body context has received little at- 
tention. Although trace formulas are applicable to inter- 
acting many-body systems, more effort has gone into de- 
veloping semiclassical descriptions of single-particle dy- 
namics in an appropriate mean field. One impressive ex- 
ception is the application of the Gutzwiller trace formula 
to the study of two-electron atoms and related three- 
body systems [^3[ Q . The main difficulty of applying the 
theory is that periodic orbits must be found for the inter- 
acting many-body system. One approach to this problem 
has been proposed in Ref. (2^] which develops a particle 
number expansion of the trace formula. 

In this paper, we consider the extension of the 
Gutzwiller trace formula to systems of noninteracting 
identical particles. The semiclassical analysis of such sys- 
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terns is quite subtle. As we shall discuss below, if the par- 
ticles are noninteracting, there is a continuous symmetry. 
(Discrete symmetries in semiclassical trace formulas are 
discussed in Refs. |26|, ||^] and continuous symmetries 
in Refs. ||, 2^, [j(J.) This symmetry has a profound con- 
sequence; although an iV-particle system in d dimensions 
can be thought of as a single-particle system in Nd di- 
mensions, one cannot simply apply the Gutzwiller trace 
formula since the presence of a continuous symmetry im- 
plies the periodic orbits of the full phase space are not 
isolated, but rather occur in continuous families. 

Recently, we presented a semiclassical formalism for 
the density of states of two noninteracting identical par- 
ticles based on an asymptotic analysis of convolution in- 
tegrals that arise in the decomposition of the semiclassi- 
cal two-particle density of states pl[ . In principle, this 
approach can be generalised to more than two identical 
particles. In this paper, we pursue a different approach 
which uses the full phase space rather than the individual 
phase spaces of each particle. We show that the formal- 
ism for continuous symmetries p9| can be used to find 
many-body trace formulas. This approach recovers our 
previous results, but can also be more easily generalised 
to arbitrary particle numbers. In addition, it is concep- 
tually cleaner than the convolution method since spu- 
rious end-point contributions from convolution integrals 
do not arise and therefore need not be explained away. 
The most important difference is that the convolution 
method cannot be used when there are interactions be- 
tween the particles whereas the analysis of this paper can 
be extended to include interactions 

It is also important to understand the effect of par- 
ticle symmetry on the semiclassical structure of many- 
body trace formulas. For noninteracting identical par- 
ticles, there are coexisting discrete and continuous sym- 
metries. While Ref. f35j] considers the symmetry-reduced 
trace formula due to the discrete permutational symme- 
try, it is assumed the periodic orbits are isolated, which 
is only true if the particles are strongly interacting (al- 
though there is a brief discussion of the noninteracting 
case). We include the appropriate continuous symme- 
tries to determine the trace formulas for the bosonic and 



2 



fermionic densities of states. 

This paper is organised as follows. In section II, we 
study the case of two noninteracting identical particles. 
We first pro vide the necessary background material in 
sections II A- II C and then give the semiclassical formu- 
lation in the full phase space. Section [II considers the ex- 



tension to N identical noninteracting particles. The sym- 
metry decomposition of the A-particle density of states is 
examined in section IV . The results of a numerical study 



of the three-particle cardioid billiard are then presented 
to illustrate and test the results of the paper. We finish 
the paper with a conclusion and several appendices. 



II. TWO NONINTERACTING IDENTICAL 
PARTICLES 

A. Quantum density of states 

The quantum Hamiltonian for two identical noninter- 
acting particles, a and b, is 



H = h{z a ) + h(z b ), 



(1) 



where z a / b denote the set of operators (£ a /biPa/f>) &ud h 
is a one-particle Hamiltonian. The full Hamiltonian (Q) 
is invariant under the unitary transformation U which 
exchanges a and b. We define the single-particle energies 
and eigenstates by 



h\j)=e j \j)- 



(2) 



Then, the two-particle energies and eigenstates are Eij — 
£i + Ej and \ij) so that 



H\ij) =E lJ \ij) 



(3) 



Accordingly, the one and two-particle densities of states 
are 
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p 2 (E) = YjiE-Eij), 



id 



(4) 



and they are related by the convolution identity 

p 2 (E) = (pi* Pl )(E). (5) 

A useful result is the relation between the density of 
states and the trace of the energy Green function or re- 
solvent. We define g(E) = Tr(G(E)), where G(E) = 
1/ (E — H) is the one-sided Fourier transform of the 
quantum propagator. In terms of the resolvent, 



p(E) = lm{g(E + ie)}. 



(G) 



and this applies for either the one or two-particle density 
of states as long as we use the appropriate resolvent on 
the right-hand side. In the limit e — > + , the exact den- 
sity of states is recovered ||. Henceforth, the ie will be 
implicit. 



B. Symmetry decomposition 

The most interesting aspect of the existence of identical 
particles is the fact that only certain states are occupied, 
the fully symmetric ones if the particles are bosons or 
the fully antisymmetric ones if the particles are fermions. 
It is important to understand how the above discussion 
decomposes when we consider the separate densities of 
symmetric and antisymmetric states. Although not ab- 
solutely necessary for the present discussion, it will be 
useful for later to introduce projection operators. As 
mentioned, the Hamiltonian (fil) is invariant under ex- 
change of the particles a and 6, an operation we denote 
by <; (leaving the particles unchanged we denote by l). 
There is a two-element discrete group comprised of these 
operations and the representation of these group elements 
in the Hilbert space (i.e. the quantum operators that 
exchange the particles) are U and I. Both of these op- 
erators commute with H. This is a simple group with 
two irreducible representations (irreps) which we iden- 
tify as the bosonic (symmetric) representation and the 
fermionic (antisymmetric) representation. Given an ar- 
bitrary state with components belonging to both irreps, 
we can project out the portion belonging to each irrep 
through the use of the projection operators |34[ 



(7) 



where the ± refer to the bosonic and fermionic irreps, 
respectively. 

In terms of these projection operators, the bosonic and 
fermionic densities of states are given as 



P± (E)=Tr(P ± S(E-H) 



(8) 



The sum of the bosonic and fermionic densities is the 
complete two-particle density of states, P2(E). The dif- 
ference is given by Ty(U6(E — H)) and expressing the 
trace in the energy eigenbasis, 

p+(E)-p_(E) = J2(ij\US(E - H)\ij) 
= J2(ji\ij)S(E - E^) 
= ^(£-2^), (9) 



where we have used the fact that U exchanges the state 
labels in the second line and the fact that Ejj = 2ej in 
the third. The final line we recognise as pi(E/2)/2 and 
thereby conclude 



p±(E) = Up 2 (E)±^ Pl 



(10) 
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C. Review of semiclassical formulation 

It is common to decompose the semiclassical density of 
states into smooth and oscillatory components. For the 
one-particle density 



Pi c (e) = pi(e)+p 1 (e), 



(11) 



where p and p denote the smooth and oscillating com- 
ponents, respectively. There is an extensive literature on 
this decomposition || . We adopt the point of view that 
one can simply use the leading-order contribution of each 
component. We do not consider the subtle issues related 
to the asymptotic nature of this decomposition (see, for 
example, Refs. 
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36[ |37|). For analytic potentials in n 
dimensions, the leading-order term for the smooth den- 
sity of states is 



Pi(e) 



1 



dzS(e — h(z)), 



(12) 



where z collectively denotes the 2n classical phase space 
coordinates and h(z) is the classical Hamiltonian (for an 
exception to this general result, see Refs. |p8|). There 
are corrections to ( jl2| ) involving derivatives of the delta 
function in the integrand. The first correction is 0(h 2 ). 
For a two-dimensional billiard, the analogous expression 
is 



, . aA y/aC 



O(e), 



(13) 



where a = 2m/ "h , A and C refer to the area and perime- 
ter, respectively and the ± refers to Neumann and Dirich- 
let boundary conditions, respectively. The third term 



K = 



1 

12tt 



&Ik{1) 



L v — 



24tt 



(14) 



is the average curvature integrated along the boundary 
with corrections due to corners with angles 0j. It does not 
actually contribute to the density of states, but rather to 
its first integral and this term will be used in section [v|. 
There are also corrections involving powers and deriva- 
tives of the curvature 
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for more exhaus- 



(see Refs. 

tive studies). Similar results hold for higher dimensional 
billiards (see Ref. M). 

The oscillating component can be written as |jj 



pi{e) « Im 

7T 



expj 



SJe) 




(15) 



where 7 labels the periodic orbits of the system, 5 7 is the 
classical action integral along the orbit and <r 7 is a topo- 
logical index ||| counting the caustics in phase space 
encountered by the orbit. A 1 is the amplitude of the 
periodic orbit and depends on the specific nature of the 
system, for example, whether the orbit is isolated or not 



and if so on its stability. For the case of isolated periodic 
orbits, 



A 1 {e) 



det(M 7 - I) 



(16) 



where T® is the primitive period of the orbit and M 7 is 
the 2(n — 1) x 2(n— 1) symplectic stability matrix on 
any Poincare section to which the orbit is transverse. Its 
eigenvalues give the stability exponents of the orbit. 

As stated above, the density of states for two noninter- 
acting particles is the auto convolution of the one-particle 
density of states (||). Formally, the semiclassical two- 
particle density of states is the autoconvolution of jll]), 
that is, 

pf(E) = (p? * pf ) (E) = p 2 (E) + p 2 (E), (17) 
where 

p 2 (E) = (p 1 *p 1 )(E), (18a) 
p 2 (E) = 2 Go x * px) (E) + (p! * fa) (E). (18b) 

The mixed term 2 (pi * pi) (E) also belongs to the oscil- 
lating component of p| c (i?). This is because an asymp- 
totic endpoint analysis of the convolution integral results 
in an oscillatory function as shown in Ref. |31| where all 
components have been evaluated and given explicit semi- 
classical interpretations in terms of one and two-particle 
dynamics which support this decomposition. We also 
showed above how the difference between the bosonic and 
fermionic densities is given by the one-particle density of 
states. Formally, we may write (analogous to Eq. (|l(J)), 



P S £{E) = \U c {E)± l -p\ c (^ 



(19) 



The formal results ( ^7[jl9| ) can also be understood from a 
semiclassical analysis in the full phase space. This anal- 
ysis is not only more fundamental, it is necessary if one 
wants to include interparticle interactions since the parti- 
cle dynamics then become coupled and we can no longer 
make use of calculations that involve the individual one- 
particle phase spaces. In the following sections, we derive 
trace formulas for the full and symmetry-reduced densi- 
ties of states from semiclassical calculations in the full 
two-particle phase space. Since we are mainly concerned 
with the extension of the Gutzwiller theory, we focus on 
the fluctuating part of the density of states. However, 
since the smooth part is important in constructing the 
complete density of states, we have provided a discussion 
of the two-particle Thomas-Fermi term (and the asso- 
ciated symmetry decomposition based on the theory of 
Ref. |lOj) in Appendix |A|. To calculate the fluctuating 
part of the density of states, we need to find all periodic 
orbits in the full phase space at a specified energy E. 
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D. Two-particle dynamics in the full phase space 

The two identical particles, a and b, evolve inde- 
pendently in their own one-particle configuration space 
which we denote as having dimension d so that the one- 
particle phase spaces are of dimension 2d. The full two- 
particle configuration space has dimension 2d and the 
corresponding phase space is of dimension Ad. We reserve 
the symbol z to collectively denote these 4c? phase space 
coordinates and will use z = (z a ,Zb) where z a /b denote 
the 2d-dimensional one-particle phase space coordinates 
of each particle. We recall that dynamics in the full phase 
space consists of each particle evolving separately in its 
own phase space. The dynamics of z arc defined through 
one-particle dynamics by $tz = (4>tz a ,4>tZb), where <j>t 
is the flow for one particle. The (nonintcracting) two- 
particle Hamiltonian is H(z) — h(z a ) + h(zb), where 
h(z a /(,) is a one-particle Hamiltonian. 

We seek periodic orbits with phase space coordinates z' 
such that $tz' — z' for some period T. This is possible 
if the two particles are on (generally distinct) periodic 
orbits with the same period. In general, two arbitrary 
periodic orbits will have different periods. However, there 
is a parameter which we can vary, namely the way in 
which the total energy is partitioned between the two 
particles. Generally, we can find an energy E a (and Eb = 
E — E a ) such that the two periods are the same. We 
will assume henceforth that there is only one energy E a 
for which there is a solution. (This assumption can be 
relaxed at the cost of heavier notation.) There is another 
way to have a periodic orbit in the full phase space; one 
particle can evolve dynamically on a periodic orbit with 
all of the energy while the other is stationary at a fixed 
point of the potential. This is discussed later. 



1. Dynamical periodic orbits 

If both particles are on periodic orbits, we call the full 
phase space periodic orbit a dynamical periodic orbit. 
We first note that such orbits occur in continuous fam- 
ilies. To see this, imagine that a full phase space pe- 
riodic orbit consists of one particle on a periodic orbit 
7o and the other particle on a distinct periodic orbit 7& 
(see Fig. [l]) and that the energy partition is such that 
both orbits have the same period T. We have complete 
freedom in specifying which points on the respective or- 
bits we choose as initial conditions. Given that we define 
t = to be when particle b is at some specified point 
on 7b, we can vary the position of particle a on "f a . By 
changing its initial position along the orbit, we map out 
a continuous family of congruent periodic orbits. 

This can be formalised as follows. We note that in ad- 
dition to the total Hamiltonian H, there is a second con- 
stant of motion J = h(z a ) in involution with H . It gener- 
ates time translations of particle a while leaving particle 
b fixed. (In fact, J can be chosen as any linear combina- 
tion of h(z a ) and h(zb) as long as it is independent of H.) 




FIG. 1: Two periodic orbits 7 a and 76 which constitute a 
periodic orbit V of the full phase space. The full Hamilto- 
nian H(z a , Zb) generates time translations for both particles 
(as denoted by the single-particle flow 4>t acting on both par- 
ticles) while the single-particle Hamiltonian J = h(z a ) gen- 
erates time translations for particle a while leaving particle 
b fixed (as denoted by the single-particle flow (j>e acting on 
particle a only). The flows generated by H and J are $t and 
$fg, respectively. A combination of such flows (cf. Eq. (po|)) 
is shown here. 

Flows generated by J are denoted by ^>g and are mapped 
in the full phase space as follows: #gz — (<j>$z a , Zb). The 
symmetry parameter 9 is conjugate to J and has the 
dimension and interpretation of time. However, since it 
only measures the evolution of particle a, it is not time in 
the usual sense and we will follow the notation of Ref . 
in denoting the paramater by 9. A combination of flows 
in H and J is 

§t^ez = {(j>t4>eZa,<l>tZb) = (4>t+ez a , 4>tZb)- (20) 

Since ^0 and $ t commute and separately conserve both 
H and J, the surface mapped out by these flows has con- 
stant H and J (i.e. H = E and J = E a ) . Starting at 
some point on the full phase space periodic orbit, flows in 
H and J map out a two-dimensional torus. This means 
there is a I-parameter degenerate family of periodic or- 
bits (the other dimension is parameterised by time and 
is present even in the case of isolated orbits). Therefore, 
we cannot use the Gutzwiller trace formula for isolated 
orbits since it will give a spurious infinity. Due to the 
continuous family, there is one fewer stationary phase in- 
tegrals to be done in evaluating the trace so that this fam- 
ily of orbits contributes 0(1/ 'VK) more strongly than an 
isolated orbit and the calculation of its amplitude must 
be performed carefully. 

For the present, we assume that there are no sym- 
metries other than J so that all periodic orbits of the 
one-particle phase space are isolated. The flow directions 
generated by H and J are stable as are the two directions 
transverse to the constant H and J surfaces. Thus, there 
are four directions of neutral stability in phase space. 
The remaining (Ad — 4) directions decompose into sepa- 
rate subspaces of dimension (2d— 2) within each of which 



5 



there are the standard symplectic possibilities for stabil- 
ity. 

In general, the leading-order contribution of one /- 
parameter family of orbits (generated by Abelian sym- 
metry) to the resolvent is p9| 



9r(E) = 



1 



1 



1 9© 


1/2 


1 dJ 


r 



det 



1/2 



x exp i 



f Sr(E) 



(21) 



This contribution is 0(l/h^ 2 ) stronger than an isolated 
periodic orbit. As mentioned above, every constant of 
motion implies one fewer stationary phase integrals and 
therefore / fewer powers of y/K in the prefactor. For 
a similar reason, there is an additional phase factor of 
— /7r/4. The total contribution to the resolvent is a sum 
over all families of periodic orbits, the capital T indicating 
that these are indeed families and not isolated orbits as 
in the more familiar Gutzwillcr trace formula. In our 
case, the sum over T can be expressed as a double sum 
over "f a and j b indicating the periodic orbits on which 
the particles are evolving. We now describe the various 
factors in (|2l]) and explain what they are in the present 
situation for which / = 1. 

The volume term, T®V® is the integral over the flows 
generated by H and J, § r dtd9, integrated over the pe- 
riodic orbit family. The time integral gives the period of 
the family T T = T la (E a ) = E /b (E b =E-E a ) = T while 
the 9 integral gives Vr = T la (E a ) since a flow in J by that 
amount returns particle a to where it began. (Hence, the 
intial phase space coordinate is mapped back to itself un- 
der the dynamics.) However, there can be discrete sym- 
metries such that a combination of flows in H and J for 
less than T restores the intial conditions. This situation 
occurs when one or both particles is on a repetition of 
some primitive periodic orbit. To see this, suppose that 
particle a is on the n a 'th repetition of its orbit while par- 
ticle b is on the n b 'th repetition of its orbit. Then, the 
torus is partitioned into n a n b equivalent segments and 
the primitive volume term is TpVj? = TrVr /n a n b . How- 
ever, the full periods are defined through the primitive 
periods by T~ fa (E a ) = n a T® a (E a ) and similarly for parti- 
cle b. Thus, T°V£ = T° a (E a a )Tl(E b =E- E a ), which is 
the product of the primitive periods. 

Mr is a (Ad — A) x (Ad — A) matrix linearising mo- 
tion on a reduced surface of section. Specifically, it is 
the section at constant (H, J, x» a , Xuf,) where X» a / b are 
chosen so that the dynamics are transverse to the sur- 
face on which they are both constant. In our case, 
this section is simply the direct product of the normal 
Poincare surfaces of section for each of the two motions 
(where one would specify the one-particle energy and 
some fixed coordinate in each case). As a result, Mr 
has a block diagonal structure since there is no cou- 
pling between the two particle spaces. We conclude that 
det(M r -I)= det(M 7a - 7)det(M 7b - I), where M la/lb 



are the stability matrices of the separate one-particle pe- 
riodic orbits and I is the appropriately dimensioned unit 
matrix on both sides of the equality. 

The anholonomy term (90/9 J) r measures the amount 
by which orbits that are periodic in the symmetry- 
reduced dynamics fail to be periodic in the full phase 
space. Suppose we vary the value of J infinitcsmally 
while keeping the total energy fixed; in our case this 
amounts to a slight change of the energy partition be- 
tween the two particles. The periodic orbit is launched 
as before with the same initial conditions except for py 
(the momentum conjugate to xn) which must be changed 
appropriately to effect the change in J. After the origi- 
nal period T, an initial phase space coordinate will not be 
mapped back to where it began, but rather infinitcsmally 
close to this initial condition. A flow in H for some extra 
amount of time At and a flow in J by an extra amount 
A9 (or vice- versa since the flows commute) closes the or- 
bit in the full phase space. The factor d<d/dJ is simply 
the ratio AO/ A J (in the limit that A J — > 0). (6 is cap- 
italised to stress that J and 9 can also be used as labels 
of families of surfaces, in which case, this factor can be 
interpreted as a Jacobian for a change of label from J to 
9.) Recall that the value of J = h(z a ) is just the energy 
of particle a. If J a — > J a + AJ a , then E b — > E b — AJ a , 
since the total energy is fixed. j a now has a perturbed 
period, T + AT 7a = T + T^AJ a while j b now has a per- 
turbed period T + AT 7b =f-T! rb AJ a , where the primes 
denote differentiation with respect to energy: 



rpl = dT -(. 

la 



de 



T' = — 

76 



dT lb (E-e) 



de 



(22) 

Let z' — (z' al z' b ) and z denote the initial and final phase 
space coordinates, respectively. Then, after the original 
period T, 



$tz' = z = (6- 



X>^-AT 7 . z' b ). 



(23) 



We need to find (At, A9) which map z back to z' . Using 
Eq. (|20|), the condition for a periodic orbit &At^Aez = z' 
implies At = AT lb and A9 = AT- 



AT 7b so that 



de =T , T , 



(24) 



The action S r (E) = S Ja (E a ) + S lb (E b = E - E a ) is 
the action of the periodic orbits in the family (all or- 
bits in T have the same action because of symmetry). 
Finally, we discuss the phase indices, /ir is determined 
from the dynamics in the symmetry-reduced surface of 
section in the same way as for isolated orbits in the usual 
Gutzwiller trace formula and following the same logic as 
above, fir = &~ la + <^-y b ■ <5r is defined as the number of 
positive eigenvalues of (dQ/dJ) T In this case, the 

anholonomy term is simply a scalar and therefore <5r = 1 
if the Jacobian is positive and Sr = if the Jacobian 
is negative. We conclude that the contribution to the 
resolvent from one family of dynamical orbits is 
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9r(E) 



if i 



'la 2 



ih 



det 



(M la -l) 



1/2 



exp i 



' S th ( E b) 



'lb 2 



expz(<5 r § -f) 



itt 



dot 



(m 76 - /) 



1/2 



Tt 



TL 



(25) 



As mentioned above, we assumed that there is only one 
energy partition such that both particles have the same 
period. This will be the case when the period is a mono- 
tonic function of energy, which is a typical situation. If 
the period is a more complicated function of energy, there 
may be further solutions and if so then one must have a 
sum over (7 , 7&) for each possible solution of this condi- 
tion, but we suppress this possibility for notational sim- 
plicity. 

We obtained Eq. ( p5| ) in Ref. j3l[] by doing a stationary 
phase analysis of the direct auto convolution of Eq. (15). 
(The phase index v in Eq. (18) of Ref. |u| has a differ- 
ent definition than 5r in Eq. (pa), but the overall phase 
is consistent in the two formulas.) The condition of sta- 
tionary phase immediately implied that the energy must 
be partitioned so that the periods of the two orbits are 
the same. The stationary phase integral then introduces 
a factor of \fh as well as the sum of the second deriva- 
tives of the actions with respect to energy evaluated at 
the stationary phase energy. This is precisely the first 
derivatives of the periods with energy. Thus, we have 
shown how these two different approaches yield consis- 
tent results. 

Either particle can execute any number of repetitions 
of its primitive orbit. If particle a executes L (a repetitions 
and particle b executes L (b repetitions, then the energy 
must be partitioned so that l la T la (E a ) = ^ 7b T 7i) (£?;,). As 
discussed above, the volume term remains unchanged as 
the square of the primitive periods due to the discrete 
symmetry in the family of orbits. The action of this orbit 
is l la S la (E a ) + l~ fh S lb (Eb), and similarly for the phase in- 
dex. The separate monodromy matrices are raised to the 
appropriate power. The anholonomy term follows from 
the discussion above Eq. (§|) as l la T^ a (E a ) + l lb T^ b (E b ). 
Assuming that the periods are monotonic in the energy, 
the index dr is unchanged. Otherwise, it depends on the 
sign of l~/ a T' + l~/ b T'- Therefore, the previous expres- 
sion continues to apply, but the actions, phase indices 
and periods are multiplied by the appropriate value of 
l la / lb and the separate monodromy matrices are raised 
to the appropriate power. This dependence can also be 
understood to be already contained in the definition of 
the various orbit properties. 

We observe that the amplitude of ( p5f ) is proportional 
to the product of the amplitudes for the single-particle 
dynamics. The trace formula for two nonintcracting par- 
ticles contains an additional prefactor of ih/VZnh, a fac- 
tor involving the derivatives of the periods with respect 
to energy (and the associated phase index S) and an addi- 
tional phase factor of 7r/4. This result generalises to cases 



where the amplitudes are not given by (|l6|). We simply 
replace the single-particle amplitudes in large brackets by 
the equivalent ones for the system under consideration. 
This can be understood by noting that the only coupling 
between the particles is as we have described and any fur- 
ther symmetry can be handled within the single-particle 
phase spaces. This conclusion can also be understood 
in the convolution picture by simply using the appropri- 
ate single-particle amplitudes when doing the stationary 
phase analysis fl31fl. 



2. Symmetry decomposition: dynamical pseudo-periodic 
orbits (DPPOs) 

As discussed initially by Gutzwiller |2(| and later in 
more generality by Robbins p7{ , in the presence of a 
discrete symmetry, the fluctuating density of states can 
be decomposed among the various irreducible represen- 
tations (irreps). For the two-particle case, this is simply 
the symmetric (bosonic) and antisymmetric (fcrmionic) 
cases. To evaluate the separate densities of states, one 
must calculate 9±{E) = Tr(P±G(E)) using the projec- 
tion operators in (Q). The first term of the projection 
operator results in the standard sum over dynamical pe- 
riodic orbits (|2^). There is a factor of 1/2 which indicates 
that this contribution is simply divided evenly between 
the symmetric and antisymmetric spectra. It is the sec- 
ond term of the projection operator which requires care- 
ful analysis. 

The oscillating part of Ti(UG(E)) can be expressed 
in terms of orbits on which particles begin at a point 
in phase space, evolve for some time T, are then ex- 
changed using the classical analogue of U with the net 
result that the particles are returned to their initial con- 
ditions. We call these orbits pseudo-periodic to distin- 
guish them from the (standard) dynamical periodic orbits 
discussed earlier. We first define the symplectic map- 
ping u corresponding to classical particle exchange as 
u{z a ,zi,) — {z b ,z a ). It has the property that u 2 is the 
identity mapping. The combination of time evolution for 
time t and particle exchange maps a phase space point 
z' = (z' a ,z' b ) to z = u^ t z' = (4>tz' b , 4>tz' a ) ■ To find orbits 
which are periodic under these combined operations, we 
require phase space coordinates z' and periods T such 
that z' — u<&tz'- Applying this combined operation 
twice, we find that z' = $2tz'. This is just the con- 
dition for a periodic orbit of period 2T in the full phase 
space without particle exchange. So we conclude that the 
initial coordinate z' must be on a periodic orbit of the full 
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FIG. 2: A dynamical pseudo-periodic orbit (DPPO) of the 
full (two-particle) phase space is constructed by placing two 
particles on a periodic orbit of the one-particle phase space. If 
E a = Eb and the particles are half a period out of phase, then 
after the combined operations of time evolution and particle 
exchange, the initial conditions are restored. 



phase space. However, this condition is still more restric- 
tive since the above considerations also imply that after 
time T, particle a must be where particle b began and 
vice-versa. This is only possible if the two particles are 
traversing the same periodic orbit, with the same energy 
and furthermore are exactly half a period out of phase. 
We shall call this a type-1 dynamical pseudo-periodic or- 
bit (DPPO). There is also the degenerate case where both 
particles begin and evolve together. This shall be called 
a type-0 DPPO and is discussed below. 

Therefore, the set of possible pseudo-periodic orbits 
is much more restricted than the set of standard peri- 
odic orbits since we have only contributions when both 
particles are executing the same dynamics. Furthermore, 
these orbits are isolated and do not come in a 1-parameter 
family. The existence of families for the standard peri- 
odic orbits is due to the freedom in specifying the rela- 
tive phases of the two motions. We no longer have this 
freedom. This immediately implies that contributions 
from pseudo-periodic orbits will be weaker by yh be- 
cause there is one more stationary phase integral to do 
than for the standard periodic orbits. (This can also be 
understood from the fact that particle exchange does not 
conserve the separate energies and so does not commute 
with J.) Therefore, the usual Gutzwillcr trace formula 
applies and we use it to determine the actions, periods 
and stabilities of these isolated pseudo-periodic orbits. 

Consider an arbitrary periodic orbit 7 of the one- 
particle phase space with period T 7 and choose some 
arbitrary initial condition on it which we shall call z' a . 
To have a pseudo-periodic orbit in the full phase space, 
we can begin at z' = (z' a , z' b — 4>T- 1 /2 z 'a)- If we now f° r 
a time T 7 /2 and then exchange the particles, we map 
z' onto itself (see Fig. |^). Therefore, the set of pseudo- 
periodic orbits in the full two-particle phase space is one- 
to-one with the set of standard periodic orbits in the one- 
particle phase space. The periods of the pseudo-periodic 
orbits in the full phase space are one-half of the peri- 
ods of the corresponding standard periodic orbits in the 
one-particle phase space. Nevertheless, when evaluating 
the trace integral we must integrate over all initial con- 
ditions on the orbit and this gives a full factor of T 7 in 
the amplitude. The actions and phase indices for the 



pseudo orbit are the same as for the standard orbit; al- 
though we integrate for only half the time, both particles 
are in motion and between them, they execute one full 
motion of the periodic orbit. The stability matrix in the 
full phase space requires careful analysis. Let M 7 be the 
stability matrix of the full periodic orbit 7 of the one- 
particle phase space and My be the stability matrix of 
the pseudo-periodic orbit 7' in the full phase space. It is 
shown in Appendix [j| that det (My —I) =4 det(M 7 - 1), 
where on each side of the equation / is understood to be 
an appropriately dimensioned unit matrix. We conclude 
that the contribution of this orbit to the oscillating part 
of Tr(UG(E)) is 



1 



rpO 

> 



2ih 



det(M 7 - I) 



: expz 



where all classical quantities are evaluated at the single- 
particle energy E/2, (Recall the symbol c denotes the 
group element that exchanges the particles.) Apart from 
the energy dependence and the factor of two in the de- 
nominator, this contribution is the same as the one from 
the corresponding primitive orbit for the single-particle 
density of states [Eqs. ( p^| ) and (|jq)]. 

As mentioned above, there is also the situation where 
both particles start at the same point on the orbit and 
evolve together. Interchanging them at the end trivially 
returns them to the same coordinates. This pseudo-orbit 
has action 25*, but should not be confused with the stan- 
dard dynamical orbit where the two particles start at 
independent points on the orbit and therefore occur in 
a 1-parameter family. The fact that we interchange the 
particles at the end ensures that the pseudo-orbit de- 
scribed here is isolated and does not occur in a family. 
The two types of orbits share the same action, but the 
standard orbit has a larger amplitude due to the differ- 
ent h prefactor and will tend to dominate. This situation 
of coexisting contributions with the same action is anal- 
ogous to a potential system with a reflection symmetry 
where there is a boundary orbit which contributes to both 
the identity term in the density of states and also to the 
reflection term. The only difference here is that the two 
types of orbits contribute with different powers of h. 

The analysis of the contribution of the type-0 DPPO 
is similar to above. We state without proof that its am- 
plitude is simply the same as the double repetition of the 
orbit, again divided by an overall factor of two as dis- 
cussed in Appendix |b| This pattern continues for higher 
repetitions, where for odd multiples of the action the par- 
ticles start T 7 /2 out of phase while for even multiples 
they start in phase and interfere with stronger (in an Ti 
sense) contributions from the standard dynamical orbits. 
Apart from the energy dependence and the factor of two 
in the denominator, the sum over repetitions is the same 
as for the single-particle density of states. 

Thus, we see that the contribution of the pseudo- 
periodic orbits to the bosonic and fermionic densities of 
states is precisely the same as the fluctuating density 
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of states of the one-particle spectrum except that it is 
to be evaluated at half the total energy (since the total 
energy is partitioned equally between the two particles) 
and should also be divided by an overall factor of two. 
In conclusion, 



P ± {E) = \U(E)± 1 -p l (f 



(27) 



consistent with (|l0|). 



E. One-particle dynamics in the full phase space 

We now discuss contributions to the resolvent from pe- 
riodic orbits in the full phase space where one particle 
executes dynamics while the other particle remains sta- 
tionary. In particular, suppose that particle a is station- 
ary at some point in phase space while particle b evolves 
dynamically on a periodic orbit. We call this a hetero- 
geneous periodic orbit. The structure of such orbits is 
qualitatively different for potential systems and billiards. 

For analytic potentials, the stationary particle must 
be at some extremum of the potential with zero momen- 
tum. In this case, the full heterogeneous orbit is isolated 
in phase space since a flow in J = h(z a ) does not map 
an initial condition z' to any new phase space point z. 
Therefore, we can use the Gutzwiller trace formula for 
isolated orbits. In billiards, the stationary particle has 
zero momentum, but it can be anywhere in the billiard. 
So rather than being isolated, the heterogeneous orbits 
appear as c£-dimensional families. This means that we 
can use the formalism of Ref. |^9| to calculate the ampli- 
tude of these orbits. 

The symmetry decomposition for heterogeneous orbits 
is trivial. Since the two particles are executing com- 
pletely different dynamics, the combination of time evo- 
lution and particle exchange, as above, can never return 
the particles to their initial conditions. This requires an 
equivalence of the two motions. Thus, the contribution 
from heterogeneous orbits is simply divided evenly be- 
tween the symmetric and antisymmetric representations 
and belongs to the pi (E) term of Eq. (|2^) . 



the full heterogeneous orbit and M a is the 2c? x 2d mon- 
odromy matrix of particle a. Since the dynamics of par- 
ticle a are locally harmonic, we can use the result for 
a (i-dimensional harmonic oscillator (see Appendix |^), 
y/\det(M a - I)\ = n j d = i2sin(w i T 7 (i;)/2). The phase 
index of this motion is simply d, one for each transverse 
harmonic degree of freedom. Thus, the contribution of 
one heterogeneous orbit to the resolvent is 



~9r(E) 



T°JE) 



det ( My 



x exp % 



SJE) 



7T 7T 
0~-y — — d — 



(28) 



where we have retained the symbol T to denote the full 
heterogeneous orbit and 7 to stress that this is the con- 
tribution from the situation where only one particle is 
evolving dynamically. There is also an identical contri- 
bution from the situation where particle b is fixed while 
particle a evolves dynamically. As before, repetitions can 
be understood to be implicit in the definitions of the ac- 
tion, period, phase index and stability matrix. 

One can also consider extrema other than potential 
minima, such as saddles or potential maxima. We can 
expand the <i-dimensional potential around an extremum 
xq as 



V(x - x ) 




(29) 



where the £ measure the deviations of x from xq. In gen- 
eral, there are d + stable directions and tf_ — d — d + un- 
stable directions. Then, the expression (|2^) is still valid, 
but the energy of the dynamically evolving particle is re- 
placed by E — V(x$), the phase factor d7r/2 replaced by 
gL|_7t/2 and the sin (wjT 7 /2) replaced by sinh (ujjT^/2) for 
the unstable directions. Finally, we note that for smooth 
potentials, the dynamical orbits give the leading-order 
contribution to P2(E) while the hetero-orbits give cor- 
rections of higher order in h. 



1. Analytic potentials 

Suppose particle b traverses a periodic orbit 7 with 
action Sy, primitive period T 7 , stability matrix M 7 and 
topological index 07 . Particle a is assumed to be station- 
ary at a potential minimum with energy E a = 0. At the 
minimum, the potential is locally harmonic with d fre- 
quencies ujj. As explained above, the full heterogeneous 
orbit is isolated and so we can simply use the Gutzwiller 
trace formula for isolated orbits. The only required in- 
formation is the monodromy matrix in the phase space of 
particle a since det(M r - I) = det(A/ a - /) det(M 7 - I), 
where Mr is the (4d — 2) x (Ad — 2) stability matrix of 



2. Billiard systems 

As mentioned above, heterogeneous orbits in a d- 
dimensional billiard occur in d-dimensional families and 
we may therefore use Eq. (|2l| ) with / = d to determine 
the appropriate trace formula. We first consider a two- 
dimensional billiard (d = 2), although the result is easily 
generalised. For this case, the orbits do not appear as 
three-tori, but rather have the topology of B x S 1 , where 
B denotes the billiard domain and S 1 is the one-torus 
associated with the dynamics of the evolving particle b 
on the periodic orbit 7. There are two constants of the 
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motion: J\ = p Xa and J2 = p Va (J = (Ji, J2)) and these 
generate flows = {x a ,y a ). Clearly, 



. f d& \ 



/ dXg dXg 



= dct 



V 



dp Va 



(30) 



since the off-diagonal elements vanish due to the fact 
that the x and y motions are uncoupled. After particle b 
has traversed the primitive orbit n 7 times, dx a /dp Xa = 
dy a /dpy a — —n~fT®(E)/m, where T®(E) is the primitive 
period of the orbit and m is the mass of the particle. (The 
minus sign indicates that a backwards flow is required to 
close the orbits in the full phase space.) This immedi- 
ately implies that the phase index 5 = 0. The stability 
matrix defined in (|2l]) in this case is simply the stabil- 
ity matrix of the motion of particle b. The volume for a 
family of such orbits is the area of the billiard and com- 
bining all of the factors, the leading-order contribution 
of a family of heterogeneous orbits T to the two-particle 
density of states is 



Pt(E) = 



S-,(E) 



'1 2 



4tt 2 



det(M 7 - I) 



(31) 



We obtained this expression in Ref. |3l| by doing a direct 
energy convolution integral of the first term of Eq. fll"3| ) 
with Eq. ( |l5|) . This once again underlines the equivalence 
of the two methods. This naturally extends to the higher 
order terms of ( |l3| ) through a more careful analysis of the 
surface corrections, but we do not pursue this analysis 
here. Also, this result generalises to d dimensions as 



d/2 



det(M 7 -7) (n 7 TO(E)) 



. d/2 



X COS 



<j~ d— 

h 7 2 4 



(32) 



where is the d-dimensional volume of the billiard. We 
stress that this is 0(l/fi d/2 ) stron gcr than an isolated 
orbit, this factor arising from the fact that this class of 
orbits appear in d-dimensional families. The contribution 
from hetero-orbits is also O^/Ti^ 1 ^ 2 ) stronger than the 
contribution from dynamical orbits. Thus, for billiards, 
hetero-orbits give the leading-order contribution to p2 {E) 
while dynamical orbits give corrections of higher order in 
H. 



III. SEVERAL NONINTERACTING 
IDENTICAL PARTICLES 

We now consider the extension to N identical parti- 
cles. The smooth term can be written as an (N — l)-fold 



convolution integral of the single-particle smooth terms 
and can also be understood as a single integral in the N- 
particle phase space. At this point, we say no more about 
the smooth term and refer the reader to the appendices 
for further discussions. For the oscillating term, there 
are again two possibilities. Either all of the particles are 
evolving dynamically, or a subset of them is stationary 
at various potential extrema (or anywhere in a billiard). 
For the first situation, the discussion closely parallels the 
two-particle case. The only nontrivial quantity to deter- 
mine is the anholonomy matrix (d&/dJ) r . We consider 
the case N = 3, but this result readily generalises. 



A. Dynamical orbits 

In an obvious extension of the notation, there are three 
single-particle phase spaces with coordinates z a , z and z c 
so that the full three-particle phase space has coordinates 
z = (z a , Zf,, z c ) and the total Hamiltonian H (z) — h(z a ) + 
h(zb) + h(z c ). Two other constants of motion which are 
in involution with H arc J a = h(z a ) and Jb = h(zb) 
and they generate time translations of particles a and 6, 
respectively, while having no effect on the other particles. 
Flows generated by H, J a and Jb are denoted by $ t , Ag a 
and VPflj, respectively. If ^ is a single-particle flow, then 
flows in the full phase space are mapped as follows: 



$ t (z a ,z ,z c ) 
Ag a (z a , z b , z c ) 
^ e b {z a , z b , z c ) 



(4>tz a , 4>tZb, (frtZc), 

(4>g a z a ,z b ,z c ), 

(z a ,4>g b z b ,z c ). 



(33) 



The periodic orbits of the full phase space (at a given 
total energy E) can be found from the one-particle peri- 
odic orbits by balancing the energy partition among the 
three particles (i.e. varying J a and Jb while holding H 
fixed) so that all the one-particle periodic orbits have the 
same period. (The result is a three-particle periodic orbit 
in the full phase space.) Imagine a slight departure from 
this equilibrium situation so that J a — > J a + A J a while 
holding Jb and H fixed. Then, 

E a -> E a + AJ a T a -> T a + AT a 

E b -> E b T b -> T b (34) 

E c -> E c -AJ a T c -> T C + AT C 

where AT a = T' a AJ a and AT C = -T' c AJ a , the primes 
denoting differentiation with respect to energy. The ini- 
tial condition z' = (z' a , z' b , z' c ) with these modified en- 
ergies (but each particle still on its periodic orbit at 
that modified energy) is not on a periodic orbit of the 
full phase space. However, it is on a generalised pe- 
riodic orbit. That is, the trajectory can be made to 
close with additional flows in (H, J ai Jb). Imagine a flow 
in H for the original period T. The orbits of parti- 
cles a and c will fail to close by the amount by which 
their period is longer (or shorter) due to the changed en- 
ergy: §tz' = (<t>-ATa,z'a,z b ,(f>-AT c z' c )- Additional flows 
in (H, J a , Jb) close the trajectory. First, a flow in H by 
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the amount AT C returns particle c to z' c : 3>at c 3>tz' = 
{4>-AT a +AT a z' a , 4>AT c z b, z ' c ). The condition for a periodic 
orbit A^e^ A8 b ®AT a ®Tz' = z' immediately implies 

A8 a = (T' a + T' c )AJ ai 

A6 b = T c 'AJ a . (35) 

We get a similar result from a deviation in (holding 
J a and H fixed) and conclude that 

fd&\_ (V a + V c T' c \ 
\d3j~\ T> Ti+Zj- W 

The determinant is T' a T' h + T' h T c + T' c T' a and is invariant 
under a permutation of the indices. Note that we could 
have chosen the two generators J a and J c and followed 
through the analogous calculation. In that calculation, 
the anholonomy matrix would be modified by permuting 
b and c in Eq. (|36|). Therefore, the eigenvalues of d&/dJ 
are not invariant. But, since the determinant is invari- 
ant, so too is the number of positive eigenvalues which 
determines the phase index S. Therefore, the final result 
is invariant. For N > 3 particles, this generalises to 




where T p is the period of the orbit on which particle p is 
residing. This can be shown by induction. 



The other factors which go into the trace formula are 
simple to determine; the discussion is very similar to the 
two-particle case and so we refrain from going into great 
detail. For N particles, flows in H and J = (Ji, Jn-i) 
map out an iV-dimensional torus. This means there are 
(N — l)-paramctcr families of periodic orbits in the full 
phase space. The total action is the sum of all the single- 
particle actions and similarly for the total phase index 
(i. The monodromy matrix is defined holding all of the 
single-particle energies constant in such a way that it is 
block diagonal among the various single-particle motions. 
The volume of the periodic orbit family is the product of 
the primitive periods. (To see this, recall that the volume 
term T-pVr = <f r dtd9id82 ■ ■ -d#Ar_i and that the primi- 
tive volume should only count distinct configurations.) 
Using Eq. @ with / = (N - 1), we conclude that the 
contribution to the resolvent from one family of dynam- 
ical periodic orbits is 



~9r(N,E) 



ih (2 7 r^)( Ar - 1 )/ 2 



p=i 



'V 2 



exp i#r? 



det(M p -J) J\T' p {E p ) 



(38) 



E 



N 



P=l T^(E P ) 



In Eq. (|3|), we have used the label p rather than the 
more cumbersome j p to refer to the periodic orbit on 
which particle p resides. We will continue to do this 
for the remainder of the section. The phase factor 8r 
is the number of positive eigenvalues of the (N — 1) x 
(N — 1) matrix (d®/dJ)r- If all of the particles are on 
distinct orbits, then there are N\ congruent but distinct 
full phase space orbits, corresponding to the choice of 
which particle to assign to which orbit. If there is more 
than one particle on the same orbit, then the number 
of combinatoric possibilities is accordingly modified. We 
take this combinatoric factor to be implicit in the sum 
over orbits and do not explicitly account for it here. 



B. Hetero-orbits 

The other possibility is that some of the particles are 
not evolving dynamically, but rather are stationary in a 
billiard or at potential extrema. Imagine that M parti- 
cles are evolving dynamically and (N—M) are fixed at ex- 



trema. Then, these heterogeneous orbits come in {M— 1)- 
fold families. In the special case where the nonevolving 
particles are stationary at potential minima, 

<7r(M, N, E) = ,9r(M, E e ) / f[ ^f} \ . 

(39) 

The evolving particles share the energy E e = E — 
EpM+i^f^li where x p denote the positions of the 
stationary particles. We recall that d is the dimension 
of the one-particle dynamics and the Wj denote the d lo- 
cal harmonic frequencies around the minimum at which 
particle p resides. As in the two-particle case, if a par- 
ticle is at a saddle or maximum, we replace the phase 
d7r/2 by d+ir/2, where d+ denotes the number of stable 
directions and replace the sin in the amplitude by sinh 
for the unstable directions. Again, there are distinct but 
congruent heterogeneous orbits in which different parti- 
cles are chosen to be on different orbits or extrema, but 
we refrain from an explicit discussion of the combinatoric 
possibilities. 
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Next suppose that (N — M) particles are station- 
ary in a (i-dimensional billiard. In addition to the 
(M — 1) independent generators that exist for the po- 
tential system, there are (N — M)d generators J q = 
(pi, . . . , Pn-m )■ The associated flows are © q — 
(qi, . . . , qjv-A/)- (Each p and q is d-dimensional.) Since 
the generators associated with the stationary particles 
also generate new orbits, the dimensionality of the or- 
bit families is / = [M - 1) + (N - M)d. The vol- 
ume term Tp V® = § r dtdOi ■ ■ ■ d6' A/ _ 1 dq 1 • • • dq w _ M = 

I? (.Ei) • ■■Tl I {E M )tif~ M) . The phase index 5 T is the 
number of positive eigenvalues of the / x / matrix 
(d&/d3)r which has a block-diagonal structure; one 
block is the anholonomy associated with the evolving par- 
ticles analogous to Eq. ( |36| ) and one block is the anholon- 
omy associated with the stationary particles analogous to 
Eq. (|30|). Thus, the contribution to the resolvent from a 
family of billiard hetero-orbits is 

${M,N,E)=3frM,E)\ n } . 

I p=M+l m / J 

(40) 

In Eqs. T is the global period (recall that the en- 

ergies of all the dynamically evolving particles have been 
partitioned so that all of the periodic orbits have a com- 
mon period) and 5r = if M = 1. As in the two-particle 
case, hetero-orbits are more important in billiards than 
in smooth potentials. Their leading-order contribution to 
p N (E) is 0(l/h {N - MKd - 1)/2 ) stronger for billiards and 
0(h^ N ^ M '' 2 ) weaker for potentials than the correspond- 
ing contribution from the dynamical orbits. 

We now make some final comments. The above ex- 
pressions apply for any of the particles executing multi- 
ple repetitions of its primitive orbit provided the energy 
is partitioned among the dynamically evolving particles 
so that all single-particle periodic orbits have a common 
period. Then, the various orbit properties which appear 
in the formulas are understood to be those for the re- 
peated orbit. The formulas written above only account 
for the contribution of a single family of orbits. The 
oscillatory part of the resolvent is a sum over all fami- 
lies: g{E) = £ r <?r(£). We mention that Eqs. (f§||c]) 
can also be obtained from convolution integrals by doing 
a stationary phase analysis of the TV-particle dynami- 
cal term and taking appropriate combinations of saddle- 
point and end-point contributions from the various cross- 
term integrals. However, the approach outlined above is 
more illuminating since it reveals the underlying struc- 
ture of the periodic orbit families. The many-particle 
trace formulas involve only properties of periodic orbits 
of the one-particle phase space. Thus, after studying a 
one-particle system, one can immediately work out the 
details of the many-particle system. This parallels the 
situation in quantum mechanics where the problem of N 
noninteracting particles in a potential is a simple exten- 
sion of the one-particle problem. 



IV. SYMMETRY DECOMPOSITION OF THE 
TV-PARTICLE DENSITY OF STATES 

If the system consists of N identical particles, it 
is invariant under Sn, the permutation group of N 
identical particles. This group has many different 
irreps for N > 2, but we only consider the one- 
dimensional bosonic/fermionic irreps which are fully 
symmetric/antisymmetric under particle exchange. We 
first introduce the projection operators j34j 

p ± = jvT^ (±1)n ^' (41) 

T 

where ± refer to the bosonic/fermionic irreps, respec- 
tively. The sum is over the group elements r of Sn which 
denote particular permutations of the particles, U T is the 
representation of the group element in the Hilbert space 
(i.e. the quantum operator which exchanges the parti- 
cles), n T is the number of 2-particle exchanges required 
to obtain r and the factor (±1) Ut is a group character. 
For fermions, the sign of the character depends on the 
number of times two particles must be interchanged. As 
before, we need to evaluate g±(E) = Tr(P±G(E)) and 
therefore Tr(U T G(E)) for each r. This is a class func- 
tion, only depending on the cyclic structure of r. 

Consider a permutation and break it up into cycles 
. For N particles, r can be decomposed uniquely into 
mutually commuting cycles; in each of these cycles, a 
subset of the particles is being permuted. An n-cycle is 
a permutation in which only n of the particles are being 
permuted. In particular, a 1-cycle corresponds to an in- 
dividual particle being left alone, a 2-cycle corresponds to 
two particles being exchanged with each other and so on. 
A general permutation r may consist of cycles of various 
sizes and also may have several cycles of the same size. 
In general, for a given r, there are v\ 1-cycles, v^, 2-cycles 
and so on. Then, the cycle structure of a class of permu- 
tations can be given as a set of integers (y\, 1/2, ... , i>n)- 
This set v labels the conjugacy classes. Two permuta- 
tions with the same v belong to the same class and thus 
have the same value of Tr(U T G). The analysis of the pre- 
vious section can be understood as being the special case 
of the identity element. To decompose the full density 
of states, one needs to determine both the smooth and 
oscillating contributions to Ti(U T G). The smooth con- 
tribution is discussed in Appendix |y. In this section, we 
examine the oscillating contribution. 

A. Dynamical cycles 

Consider first the case for which all particles are evolv- 
ing dynamically. A group element t consists of m T cycles, 
a given cycle k consisting of interchanging nk particles. 
As in the two-particle case, particle interchange does not 
commute with all of the single-particle energies and so we 
do not expect periodic orbit families of dimension (TV— 1). 
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FIG. 3: A specific permutation of 6 particles is decomposed 
into three dynamical cycles. Each of the particles belonging 
to a particular cycle are on a periodic orbit of the one-particle 
phase space with T x {E a = E b = E c = £i)/3 = T 2 {E d = E e = 
E 2 )/2 = T 3 (E f = E 3 ) = T. 



However, for each cycle, there is a generator Jk which is 
the sum of the single-particle Hamiltonians of the par- 
ticles involved in that cycle and is preserved under the 
action of the group element r. These generators com- 
mute with each other and with the total Hamiltonian H . 
However, this is not an independent set since ^2 k Jk = H. 
There are (m r — 1) independent commuting generators 
other than the full Hamiltonian and so we expect peri- 
odic orbit families of this dimensionality contributing to 
Tr(U T G). 

We seek structures in the full phase space which are 
invariant under the combined operations of time evolu- 
tion (for time T) generated by H and particle exchange 
as specified by r. Clearly, this is only possible if all par- 
ticles of a given cycle k are on the same periodic orbit, 
7fe. They must all have the same energy, which we shall 
call Ek and then Jk = nkEk- For example, imagine that 
particles a, b and c constitute a 3-cycle. Starting with 
particle a at some arbitrary point on a periodic orbit 7 
of the one-particle phase space, particle b an amount T 7 /3 
ahead of it and particle c an amount T 7 /3 behind. Then, 
after a time T = T 7 /3, a — > b, b — > c and c — > a. How- 
ever, the group element r = (acb) maps a — > c, c — > b 
and b — > a simply undoes this change and the original 
configuration is restored. Such a cycle is shown at the 
left of Fig. H We then imagine that for every cycle com- 
prising t, there is a train of particles with identical ener- 
gies traversing a periodic orbit of the one-particle phase 
space. Each particle completes (l/nfc)th of the full mo- 
tion on the periodic orbit. 

The logic then is as follows. We assign each cycle a 
periodic orbit 7^. (We will henceforth label the orbit 
properties using the subscript k rather than the more 



cumbersome 7^.) We partition the energy (i.e. the val- 
ues of Jk) so that the values of the periods Tk/nk are all 
the same; this quantity we denote by T. After time T 
and permutation r, the resulting structure is guaranteed 
to be globally periodic in the full phase space. Such an 
orbit comes in an (m r — 1) degenerate family which can 
be understood as follows. For each cycle, it is enough to 
specify the initial condition of one particle after which we 
know the initial conditions of all the other particles. We 
choose the initial condition of the first particle arbitrarily 
for the first cycle. The first particle of the other (m T — 1) 
cycles can then begin anywhere on their respective orbits 
(this constituting the dimensionality of the family). We 
also understand this from the fact that starting at the 
arbitrary initial condition, flows generated by any of the 
(m T — 1) generators Jk map out a surface of this dimen- 
sionality. Together with a flow in H, the periodic orbit 
surface is a torus of dimension m T . 

For the symmetry decomposition (involving the dy- 
namical orbits) of a two-particle system, it was noted 
that there were contributions from higher multiples. For 
instance, one could start both particles at the same point 
on a single-particle orbit, let them evolve for a full period 
and then interchange them. There is an analogous struc- 
ture in the A^-particle case. We can allow the particles to 
execute a fraction h/rik of an orbit as depicted in Fig. |^. 
As before, the additional factor Ik can be absorbed into 
the definitions of the various classical parameters. 

The contribution of an m r -torus of orbits can now be 
inferred from our previous work. The only detail is in the 
determination of (d®/dJ). It is as in (37), but with the 
understanding that the sum/product over orbits should 
be replaced by a sum/product over cycles. These become 
equivalent in the identity contribution which was consid- 
ered there. Also, since the anholonomy term measures 
deviations away from global periodicity arising from a 
change in the energy partition (now among the cycles), 
T' should be replaced by T' k /n\. A factor of 1/nfe comes 
from the fact that the energy of the cycle must be divided 
evenly among the nk particles belonging to that cycle. A 
second factor of 1/ n% comes from the fact that the orbit 
has time Tk/nk for the anholonomy to evolve. (Note that 
if this orbit is a multiple repeat, then it is understood that 
T'k = hTk > where T k is the primitive period.) The en- 
tire contribution should also be divided by \\ k rife arising 
from the monodromy matrix as discussed in Appendix 
This last fact is the generalisation of the factor of \J2 ap- 
pearing as a prefactor in the second term of Eq. ( |27| ) for 
the two-particle case. Therefore, the contribution from a 
family of dynamical cycles to Tt(U t G) can be written as 
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FIG. 4: (Left) The same type-1 dynamical 3-cycle of Fig. ||. (Middle) A type-2 dynamical 3-cycle. The same periodic orbit, 
but each particle executes two-thirds of the complete motion and the net action is 2S. (Right) A type-0 dynamical 3-cycle. 
Each particle executes one complete motion and the net action is 3S. 



in (2iriti)( r - 



-l)/2 



n 



T°(£ fe )exp*( 



S k (E k ) 



exp i5% 



det(M fe - 7) ^m(E k )\ 



(42) 



where Mk is the stability matrix for a full cycle k (cf. 
Appendix |b|) . Notice that the contribution of the group 
element for which all of the particles belong to the same 
cycle is proportional to pi(E/N) |fts|| . 



B. Hetero-cycles 

It is also possible that Tr (U T G) has a contribution from 
cycles where some particles are fixed (either at extrema 
of the potential or anywhere in a billiard) while others 
are evolving dynamically. Let s denote the number of 
cycles that are stationary and e the number of cycles that 
are evolving dynamically. Then, s + e = m T . To have 
such a contribution to the oscillating component, group 
elements must consist of two or more cycles since those 
that consist of only one n^-cycle will either contribute to 
Eq. if they are dynamical cycles or to the smooth 
part (cf. Appendix ^) if they are stationary cycles. In 
addition, we require at least one cycle to involve particles 
which are evolving dynamically (e > 1) and at least one 
cycle to involve particles which are stationary (s > 1). 
Thus, hetero-cycles are cycles for which 1 < e < N and 
1 < s < N. 

For potentials, the dimension of a family of orbits is 
then (e — 1) since only the generators associated with dy- 
namical cycles generate new orbits. The stationary cycles 
simply contribute their monodromy matrices and phase 
indices and otherwise play no essential role. Eq. 
holds for the particles which are evolving dynamically, 



but m T is replaced by e and the energy associated with 
the e dynamical cycles, E e , is the total energy minus the 
sum of the potential energies of the stationary particles. 
For a potential minimum, the contribution of one family 
of hetero-cycles to Tr(U T G) is 

n aL\ ■ 

k=e+i n j= i2sm^ k — J J 

(43) 

where the ujj k denote the local frequencies around the 
potential minimum at which the particles of cycle k re- 
side. If this cycle of particles is actually at a saddle or a 
maximum, the final factor is modified as in the discussion 
below Eq. d|). 

For billiards, the previous relation holds for the dy- 
namical cycles, but the product over stationary cycles 
is modified. As explained below, the dimension of the 
orbit families is [(e — 1) + so?] since the generators asso- 
ciated with stationary cycles also generate new orbits. If 
there are s stationary 1-cycles, these generators and their 
flows are J = (pi, . . . , p s ) and = (qi, . . . , q„), respec- 
tively. (There are sd components since each and q^ 
is ci-dimensional.) In fact, this is true regardless of the 
number of particles belonging to the stationary cycle. At 
first, this may seem incorrect since longer cycles will in- 
troduce additional generators because they involve more 
particles. However, this larger set of generators is not 
an independent set. To see this, recall that the particles 
involved in a stationary cycle can be anywhere in the bil- 
liard. If the cycle is not a 1-cycle, but rather an nfc-cycle, 
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the combined operations of time evolution and particle 
exchange will not restore the initial configuration unless 
all the particles involved in that cycle possess the same 
phase space coordinates. More formally, a stationary re- 
cycle possesses a set of generators, J = (pi, . . . , p„ fc ) and 
associated flows S — (qi, . . . , q« fc ), where each and 
has d components. However, after the specification of a 
single (p;,q;) pair, all others are uniquely determined: 



Pi = P2 



p„ fe and qi 



q2 



Thus, 



one independent set of generators is J = Pi/rik, where 
Pi is the momentum of the ith particle of the station- 
ary cycle. The factor of l/n& comes from the fact that 
the momentum of the cycle must be equally partitioned 
among the nk particles belonging to that cycle. We note 
here that it is not necessary for stationary particles of 
distinct cycles to have the same phase space coordinates. 
Thus, for a d-dimensional billiard, the contribution to 
Tr (U T G) from a family of hetero-cycles is 



g*(e,m T ,E) = gt(e,E)\ J[ 



l d exp 



k=e+l 




'2nhn k T k 



(44) 



Eqs. ( f43f]44| ) are the most general formulas of the paper. 
We allow for any amount of particle permutation and 
any number of particles can be evolving while the rest are 
stationary. Each cycle can involve an arbitrary repetition 
of the primitive motion. As before, if e = 1, then 5 = 0. 
Hetero-cycles in billiards are 0(1/V i ' 2 ) stronger than 
in smooth potentials. For potentials [billiards], hetero- 
cycles are 0(h s/2 ) weaker [0(l/h {d ^ 1)s/2 ) stronger] than 
dynamical cycles. Thus, the most significant structures 
(in an H sense) are dynamical cycles for smooth potentials 
and hetero-cycles for billiards. 



V. NUMERICAL EXAMPLE: THE 
THREE-PARTICLE CARDIOID BILLIARD 

To illustrate the use of the trace formulas derived 
above, we study a system of three noninteracting identi- 
cal particles in a two-dimensional cardioid billiard. In a 
billiard, classical orbits possess simple scaling properties. 
For instance, the action and period of an orbit 7 are 



S 7 (e) = V2meL 



TJe) = S'Je) 




(45) 



For this reason, it is natural and convenient to analyse 
the length spectrum of the various trace formulas. Thus, 
for our analysis, we shall compare Fourier transforms of 
quantum spectra with their semiclassical approximations 
in the reciprocal space of orbit lengths, L. In recip- 
rocal -L-space, we expect signals at the lengths of the 
full periodic orbits of the three-particle system. In the 
subsequent analysis, peaks in the various length spectra 
are identified with particular periodic orbits of the full 



classical phase space. We first consider the total (three- 
particle) density of states for the cardioid system and 
then study its decomposition among the irreps of 53. 



A. Total density of states 

1. Quantum mechanics 

The analogue of Eqs. (|J) and (Eh for the quantum 
three-particle density of states is 

p 3 (E) = ^8(E-(ei + e s + e k )) 

= p 1 {E)*p 1 (E)*p 1 {E). (46) 

In fact, this relation applies even if the particles are not 
identical where the full density is still the convolution of 
the three distinct single-particle densities. We construct 
the three-particle spectrum by adding the energies of the 
one-particle spectrum. (The billiard has a reflection sym- 
metry which implies that all the single-particle states are 
either even or odd; this symmetry should not be con- 
fused with the symmetry due to particle exchange.) In 
the subsequent analysis, we work exclusively with the 
odd-parity one-particle spectrum. We include the first 
500 single-particle energies which allows us to construct 
the first 19,317,062 energy levels representing all three- 
particle energies less than 2.8148 x 10 3 . (The spectrum 
is truncated at -E max = 2ei + €500 to ensure there are no 
missing levels.) It is possible to improve the resolution in 
L-space by truncating the spectrum at a higher energy. 
But, this would require a precise spectrum since there 
is a rapid increase in the number of three-particle levels 
with energy and errors accumulate. 



2. Weyl term 

Using the identity contribution from Appendix the 
smooth three-particle density of states is just the 2-fold 
convolution integral of the smooth single-particle density 
of states: 



p 3 {E) = p 1 {E)*p 1 (E)*p 1 {E). 



(47) 



For a two-dimensional billiard, we use Eq. ( |l3| ) for pi(E). 
After performing the necessary integrations (ignoring 
terms 0(l/h 3 )), the three-particle smooth term is found 
to be 



a 3 A 3 

1287T 3 

3 



E 



a 5 / 2 A 2 C 



£3/2 



—a 

2 



327T 3 

AC 2 A 2 JC 



128tt 2 16ti- 



(48) 



E. 



For the odd-parity single-particle spectrum of the car- 
dioid, A — 37r/4, C = 6, and JC — 3/16. Some of the 
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contributions of the higher-order terms of p^{E) can be 
calculated, but it is formally meaningless to include them 
since there are corrections of the same relativ e ord er in 
h that are not known. The terms that are 0(v a 3 E) and 
0(aE°) can be computed numerically. 



This leads to a simple linear system which can be solved 
to give 



E; 



T2 

L 1i 



LI + LI 



E 



(52) 



3. Hetero-orbits 

For three particles in a two-dimensional billiard, there 
are two types of heterogeneous orbits. The first type oc- 
curs when one particle is on a periodic orbit while the 
other two particles are stationary. These orbits come 
in 4-parameter families. The trace formula is obtained 
by using Eq. @ with M = 1,N = 3. For the situa- 
tion where particles a and b are stationary and particle c 
evolves on the orbit 7, the leading-order contribution to 
h{E) is 



pf{E) 



a 3/2^2ffl/2 
8^ 



E 

7 



det(M 7 - 1) 



The second type of hetero-orbit arises from the situa- 
tion where only one particle is stationary while the other 
two evolve on periodic orbits. For instance, particle c is 
stationary while particle a evolves on j a and particle b 
evolves on 7^. Using formula (|l0|) with M = 2, N = 3, 
we conclude the leading-order contribution to p^ (E) from 
these hetero-orbits is 



™5/4/|pl/4 [ 

#w=^e n 



(2rrr' 
(Ll +L 



L 1p 



,.-y 6 ^P=a,6W|det(M 7p - 7) 

2 \-3/4 / / 7T 37T 



la ' lb) 



cos v aELr — cr - 



(50) 



+ cr- 



where L T = L^ a + L? /b and a r = {a 7a , ^ lbl 

For the total density of states, both formulas are mul- 
tiplied by a factor of 3 since there are three identical 
contributions depending on the choice of which particle 
is evolving and which is stationary. Higher order contri- 
butions can be obtained using the convolution formalism 
and the results are given in Appendix 



for i = a,b,c. We can now proceed to compute each of 
the quantities involved in formula (|38|). The anholonomy 
term (cf. Eq. (§€ 



dot 



d® 



rrif rpf ■ rjif rrif . rpf rji! 

la lb lb 7c 7c la 



16 £2 « L 2 E 3 ■ 

la lb lc 



(53) 



In addition, Tr < and this implies the phase fac- 

tor dr = 0. Then, the three-particle dynamical term can 
be written as 



(49) Pi{E) 



X cos 



e if n 



\ 



p—a,b,c 



det(Af 7p - 1) 



( 1 n n 

[VaELr-ar--- 



(54) 



where Lr = 

(°7o "I" °7b + °"7e) - 



Ll and or 



5. Numerics 



We first mention that for billiards, it is common to ex- 
press the density of states in terms of the wavenumber k, 
where e = k 2 /a so that p(k) = 2kp(e)/a. This is conve- 
nient since k is conjugate to the periodic orbit lengths L. 
Therefore, our numerical results will be quoted as func- 
tions of k with the understanding that these functions 
have been converted to the k domain from the energy 
domain using the Jacobian relation above. This will al- 
ways be the case when the argument is k. As well, for all 
numerical comparisons, a (and Ti) are set to unity. 

We compare the Fourier transform of the oscillatory 
part of the density of states 

F| c (i) = Hh(k)} 

= T{3pf(k) + 3pf(k) + p d 3 (k)} (55) 



4- Dynamical orbits 



with its quantum mechanical analogue which we define 
to be 



To use formula ( 38 ) , we must first determine the ener- 
gies that satisfy the following conditions: 



T la {E a )=T. 



E„ 



ib 
E h 



+ E C 



T. 



lc( E c), 

E. 



(51) 



(56) 



In Eq. (p6|), the first term is the quantum three-particle 
density of states, ps(k) = J2i 8 (k — ki), where the su- 
perindex I denotes a triplet of integers (i,j,k). The 
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FIG. 5: Some of the shorter periodic orbits of the cardioid in 
the full domain. The label of each orbit includes the number 
of reflections and also a letter index to further distinguish it. 
The asterisk designates a self-dual orbit j42|. The two orbits 
*8b and *10b reflect specularly near the cusp, contrary to ap- 
pearances while the orbit 4a misses the cusp. From Ref. [|42[| . 



subtracted term is the smooth density of states as de- 
termined from Eq. (Q) . The oscillatory part has contri- 
butions from hetero-orbits (^9-50) and dynamical orbits 
(|54|). In all formulas, 7 P are periodic orbits in the fun- 




damental domain (i.e. the half-cardioid) and are 
their primitive lengths. Orbit properties are discussed 
in Refs. (42], and some of the shorter geometrical or- 
bits arc shown in Fig. ||[ The stability matrices of the 
Gutzwiller amplitudes are computed using the standard 
procedure for the stability of free-flight billiards (see, for 
example, Ref. 0). We define the Fourier transform 



Hp(k)} 



w(k)exp{ikL)p(k)dk (57) 



as a function of the conjugate variable L. Here, w(fc) 
the three-term Blackman-Harris window function fl44|] 



is 




ko < k < kf 
otherwise 

(58) 

with (a ,a 1 ,a 2 ) = (0.42323,-0.49755,0.07922). We 
choose kg and kf so that the window function goes 
smoothly to zero at the first and last eigenvalues of the 
three-particle spectrum. Numerical integration of (pq) 
and (Bq) using this integral operator is displayed in FigTlfl. 
In the semiclassical transform, a total of 212 periodic or- 
bits (including multiple repetitions) were used. 

We observe good agreement between quantum and 
semiclassical results for L < 7. In fact, it is difficult to 
distinguish between the two curves. For this reason, we 
plot the difference between them in Fig. [?]. Clearly, the 
errors are small with respect to individual peak heights. 
Furthermore, the errors are largely due to hetero-orbit 
contributions. This can be understood by considering 
the first three structures in L-space. The first structure 
(L w 2.60) is due to a type-1 hetero-orbit where two par- 
ticles are stationary and one evolves on 7 = |(*2a). The 



FIG. 6: Fourier transform of the oscillatory part of the three- 
particle density of states for L < 9. The solid line is the 
transform of the quantum three-particle spectrum ([36|) and 
the dashed-dotted line is the transform of the combined semi- 
classical three-particle trace formulas (^). Each structure is 
due to one or several periodic orbits of the full phase space. 



second structure (L 3.67) is from a type-2 hetero-orbit 
where one particle is stationary and two particles evolve 
independently on the same orbit 7 = |(*2a). The third 
structure arises from the interference between a type-1 
hetero-orbit (L w 4.62) where one of the three particles 
is on 7 = ^(*4b) and a dynamical orbit (L « 4.50) where 
all three particles evolve independently on 7 = i(*2a). 
We see that the first and third structures have similar er- 
rors and thus conclude that the error introduced from the 
dynamical term is much smaller than that from the het- 
erogeneous terms. All other L-space structures arise from 
the interference of many orbits and can be accounted for 
in a similar manner. For L > 7, the discrepancies are 
more significant and mostly due to the problematic or- 
bits 7 = 4a and 7 = |(*10b) that are not well isolated 
in phase space and pass close to the cusp of the cardioid 
(cf. Table |). These orbits have inaccurate Gutzwiller 
amplitudes for reasons explained in Refs. pR M. 



B. Symmetry decomposition 

Due to the identical nature of the particles, the eigen- 
states of H can be classified according to the irreps 
of S3, the permutation group of three identical parti- 
cles. Each group element belongs to one of three classes 
((3, 0, 0),(1, 1, 0),(0, 0, 1)) based on the cycle structure of 
that element. Thus, there are also three irreps. These 
are the symmetric (trivial) irrep A+, the antisymmetric 
irrep A- and the two-dimensional mixed-symmetry irrep 
£. (Sjv always possesses exactly two one-dimensional ir- 
reps regardless of the size of N > 1.) The character table 
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FIG. 7: Fourier transform of the difference between the quan- 
tum and semiclassical density of states for L < 7. The upper 
and lower windows show the real and imaginary parts, respec- 
tively. 



1. Quantum mechanics 

The total three-particle density of states can be de- 
composed into symmetry-reduced densities px{E), each 
belonging to an irrep X of S3 : 



p 3 (E) = p + (E) + p_(E) + p £ (E). 



(59) 



Each partial density may be obtained by projection 
Pi(E) = Tr(Pz5(E—H)), where the operator Px projects 



onto the irrep X \ 
eigenbasis as in (£ 



%. Expressing the trace in the energy 
J7 the symmetry-reduced densities are 



P±(E) 



7P1 



Ps{E) = g 



(60) 
(61) 



To understand how the cross term arises in Eq. (|6 
consider the contribution from t = (a)(bc) 6 (1, 1, 0) 

J2{ijk\U T 5(E - H)\ijk) = J2( ik ^ k ) 5 ( E ~ E »k) 

i,j,k 
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7i 
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73 


7.5637 


i(*2a) 


4a 




7.5650 


±(*2a) 


if 10b) 




7.9975 


±(*2a) 


±(*2a) 


4a 


7.9987 


±(*2a) 


|(*2a) 


i(*10b) 


8.4731 


i(*4b) 


4a 




8.4742 


§(*4b) 


if 10b) 




8.8011 


2a 


4a 




8.8022 


2a 


if 10b) 




8.8624 


±(*2a) 


§(*4b) 


4a 


8.8636 


±(*2a) 


§(*4b) 


iriob) 



TABLE I: Some of the orbits responsible for numerical dis- 
crepancies. The first column gives the length of the periodic 
orbit r in the full three-particle phase space while the other 
columns specify the constituent periodic orbits 7^ of the one- 
particle phase space. (Type-2 hetero-orbits involve only two 
orbits since one of the particles is stationary.) 



for S3 is given below. Numbers in front of class labels 
indicate the number of elements in that class. 



Y,S(E~(e l + 2e j j) 



1 



*Pi(E). (62) 



The other contributions can be found in a similar man- 
ner. We could compute each partial density separately 
for comparison with the numerics, but it is more illumi- 
nating to isolate the contribution from each symmetry 
class by inverting the above system of equations: (we ig- 
nore the identity class which reproduces the total density 
of states) 



P(i,i,o)(£) = P+ {E)-p_{E) 
1 



2 P1 



P(o,o,i)(E) 



-(E) 



l ) !•/„(£). 



■ p.(E) - -p £ (E) 



(63) 



(64) 



where p( ljl )(^) and p( ,o,i) (-^) denote the densities be- 
longing to the class of two-particle and three-particle ex- 
changes, respectively. From Eq. (|64|), we note that the 
contribution of the longest cycle is directly related to the 
one-par ticle d ensity of states as discussed at the end of 
In the following sections, 



section 



IV A 



we discuss the 
semiclassical decomposition of each partial density into 
smooth and oscillatory components: 



pf(E)=p x (E)+p x (E). 



(65) 
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1(3,0,0) 


3(1,1,0) 


2(0,0,1) 


A+ 
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A- 
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-1 
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2 





-1 



2E 2 = H = E => E 2 = E/2), the result has the structure 
of the leading-order term of hpi{E/2) * Pi(E). There is 
also the situation for which the 2-cycle is stationary and 
the 1-cycle is dynamical. Using Eq. (Q) such that k = 1 
is the dynamical cycle (just a standard periodic orbit of 
the one-particle phase space) and k = 2 is the 2-cycle, 
the result has the structure of the leading-order term of 
\ P \{Ej2) * p\{E). Group elements r S (0,0,1) consist 
of single 3-cycles. Therefore, there are no contributions 
from this class. To summarise, we have shown that 



TABLE II: Character table for S 3 . 



2. Stationary cycles 

If all particles being permuted are fixed, the cycles are 
stationary and contribute to px(E). Using the results of 
Appendix [d| it can be shown that the smooth densities 
belonging to each irrep are given by Eqs. (|6C|-|6l|), but 
with p replaced by p. Thus, 



0(1,1.0) C^) 



E 
~2 

a 2 A 2 
16tt 2 



-E 



*Pi(E) 

3 / 2 AC(l 



a 



V2) 



16tt 2 




3AJC 



E 



(66) 



P{o,o,i)(E) 



E 



aA 
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87T y/E 



3K8(E) 



.(67) 



Note that in Eq. (|66j), we have ignored terms that are 
0(l/ti) since some of the contributions at this order can- 
not be calculated exactly. These terms can be computed 
numerically, but are insignificant for our analysis. 



3. Hetero-cycles 

The leading-order cycles are the three 1-cycles of the 
identity class. If one 1-cycle is stationary and the other 
two 1-cycles are dynamical (s = 1, e = 2), then the result 
from Eq. ( |l4| ) is identical to /5p(2, 3,E). If instead two 
1-cycles are stationary and one 1-cycle is dynamical (s = 
2, e = 1), then Eq. @ reduces to p£(l, 3, E). There are 
three contributions of each type. 

The first correction is from permutations r € (1, 1,0) 
that consist of one 1-cycle and one 2-cycle. There are 
two such contributions. The first one is from hetero- 
cycles for which the 1-cycle is stationary and the 2-cycle 
is dynamical. Using formula ( fl4| ) such that k = 1 is 
the 1-cycle and k = 2 is the 2-cycle (rii = 1, n2 = 2, J 2 = 



±3 ( - Pl 



[3%\E)+3&(E) 

E\ , , 1 (E 
2 )*h(E) + ~- pi - 



Px{E) 



p h £ (E)=2[pf(E)+pf(E)}. 



(68) 



(69) 



4- Dynamical cycles 

The leading-order contribution to p± (E) comes from 
the identity element i = (a)(6)(c) which consists of 
three 1-cycles (m r = 3; J\ — h a ,J 2 = hi,,,!^ = h c ; 
J2k^ k = H). Thus, there are 2 independent commut- 
ing generators other than H and so we expect periodic 
orbit families of dimension 2. Using Eq. (fl2] ) and the 
fact that 1-cycles are equivalent to periodic orbits of the 
one-particle phase space, we find the leading-order term 
oipi(E) is p d 3 (E)/6. 

The next contribution is from permutations r g 
(1,1,0). There are three elements in this class each 
consisting of one 1-cycle and one 2-cycle (m T = 2; 
k = l,ni = 1; k = 2,n 2 = 2). Then, for r = (ab)(c), 
J i = h Cl J 2 = h a + hb and similarly for the other ele- 
ments in this class. Thus, there is only one independent 
generator (other than H ) and we expect one-dimensional 
families. Using Eq. (|42|), we find this contribution has the 
structure of a two-particle density. The 1-cycles (k = 1) 
are assigned to 71 and the 2-cycles (k = 2) to 72, where 
7x/2 are any periodic orbits of the one-particle phase 
space. Then, all cycle properties are those of the cor- 
responding orbit (cf. the 1-cycle and 2-cycle of Fig. ||; 
note the repetition l 2 /n 2 = 2/2 = 1 which denotes the 
case where the particles of the 2-cycle evolve together 
is not shown). Multiple repetitions of the 2-cycle are 
either fractions (if l 2 is odd) and correspond to type-1 
DPPOs or integers (if l 2 is even) and are type-0 DP- 
POs (cf. the classification used in section [I D 2j ). The 
generators J\ = n\E\ = E\ and J 2 = n 2 E 2 = 2E 2 are 
the energies of the particles involved in the 1-cycle and 
2-cycle, respectively (particles of the 2-cycle each have 
energy E/2). Thus, the final form is structurally equiv- 
alent to \px(E/2) * pi(E). 

The two group elements r € (0,0,1) each consist of 
one 3-cycle (m T = 1, A; = 1, n\ = 3) which implies there 
are no generators independent of H and thus the orbits 
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are isolated. As before, cycle properties can be mapped 
to those of an orbit of the one-particle phase space (cf. 
the 3-cycle shown Fig. [|; lk/nk — h/ n i = 1/3,2/3,3/3; 
higher repetitions l\jn\ — l\j3 would have action liS, 
phase index l\a, stability matrix M l1 , where S,a,M are 
the properties of the primitive orbit to which the cycle 
is assigned). The energy E\ in Eq. ( |42] ) is the energy of 
each particle involved in the 3-cycle (k = 1) and since 
H = Ji = 3E X = E, it follows that E x = E/3. Thus, the 
result has the structure of a one-particle trace formula, 
but it is evaluated at E/3 and has a cycle structure pref- 
actor of 1 /3. Including the prefactors from the projection 
operator, we conclude that 



P%(E) = - 



pUe) 



2p d 3 (E) 



E 

3 Pl {j 



(70) 
(71) 



We stress that even though the correction terms have 
structures equivalent to one- and two-particle densities, 
they are in fact contributions from the dynamical cycles 
of the full three-particle phase space. 



5. Trace formulas for the two symmetry classes 

Combining the results of Eqs. j6^-[7l|), the fluctuating 
densities for the two nontrivial symmetry classes are 

P(i,i,o)(£) = p+(E)-p-(E) = ~p 1 f!pJ*p 1 (E) 



+ 



pa,i )0 )(£) + [p(W^) + p\1 i,o) ( E ) 

$.1,0) (£)+P^i,o)(£) 



(72) 



and 



P(o,o,i)(E) = p + (E) + p_(E) - -p £ (E) 

\(E 

= 3 Pl 3 



(73) 



The leading-order term of 1 ^ (E) is given by 
/5(i.i,o)(£) 



(Z°/2L 7 ) 



4tt 2 
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det ( M 1 — I 







- 7T 7T 

x cos v aEL~, — 

7 7 2 2 



(74) 
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The first term of (Q) is the contribution from two par- 
ticles being stationary at the same point in the billiard 



(i.e. a stationary 2-cycle) while the third particle evolves 
on a periodic orbit (i.e. a dynamical 1-cycle). The sec- 
ond term is the contribution from one particle being sta- 
tionary (i.e. a stationary 1-cycle) while the other two 
particles evolve on a periodic orbit (i.e. a dynamical 2- 
cycle). Higher-order contributions from hetero-cycles can 
be worked out. These are included in the numerics, but 
we do not write them out explicitly here (see Appendix 
|e|). The contribution from the dynamical cycles as de- 
termined above can be written as 

I* .„ \ 



^i,i,o) (E) 



v 3/4 



(2tt) 3 / 2 E 1 / 4 



E 



n 



V 



det M. 



(m 7< -z) 



[2(2L 7i +L 72 )]- 1/4 cos 



\l OiE 7T 7T 



(75) 



where L i2 = y2i 7i + L 72 and <ji 2 = (<r 7l + c 72 ). To 

understand how this result is obtained, recall the struc- 
ture of the dynamical cycles in this class. Each full cycle 
consists of one 1-cycle and one 2-cycle. The total energy 
is partitioned among the three particles such that the 
periods of the cycles are the same. Suppose the 1-cycle 
and 2-cycle are associated with the orbits 71 and 72, re- 
spectively. The energies E\ , E2 are determined from the 
periodicity condition 

1, 
2~ 

Using the usual relations for actions and periods in a 
billiard (Hq), one can show that 



T 11 (E 1 ) = -T 12 (E 2 ) 



(76) 



Ei = 



2LI 



2 ^ 



L 72 



E, E 2 = 



2L 71 



E, (77) 



where E\ is the energy of the particle of the 1-cycle, 2E 2 
is the total energy of the particles involved in the 2-cycle 
(each of them has energy E 2 since their energies must 
be equal) and E — E\ + 2E 2 is the total energy of the 
three-particle system. The 2-cycles are similar to the 
DPPOs of a two-particle system (cf. section II D 2) and 
we shall use the same classification scheme for all 2-cycles. 
The trace formula for p( ,o,i)(^') is a one-particle trace 
formula except that L 7 - - ! 7 ' 1 
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6. Numerics 

We first consider the class (1,1,0). We compare nu- 
merically the length spectrum of the dynamical cycles 

Ftl m {L)=T{p d (lxfi) {k)} (78) 

with its quantum analogue 

F (Tl 0)( L ) = • ?7 {P(l.l,0)(fc) ~ P(1,1,0) (*0 ~ £(1,1,0) 

(79) 
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FIG. 8: Length spectrum for the class (1,1,0). Quantum 
(solid line) and semiclassical (dashed-dotted line) results for 
L < 7. Each peak is due to a dynamical cycle of the full 
three-particle phase space. 




FIG. 9: Length spectrum of the hetero-cycles for the class 
(1,1,0). The upper and lower windows show T {pji 1 1 o)(fc)} 

and •^{/5(i ilj0 )(fc)}) respectively. 



We construct the first 241,080 levels of P(i,i,o)(fc) using 
the first 1000 single-particle energies. The smooth term 
is computed from Eq. ( |66| ) using the billiard parameters 
(as above) for the odd spectrum. Trace formulas are 
computed using geometrical orbits with length L < 10. 
The result is shown in Fig. |8[ 

We now examine some of the L-space structures of 
Fig. ||. The first peak (L « 3.18) is due to the dynami- 
cal cycle where both 1- and 2-cycles are on the primitive 



TABLE III: Some dynamical cycles of the three-particle car- 
dioid billiard for the class (1,1,0). The first column gives the 
position of the peak in L-space arising from the dynamical 
cycle while the third and fourth columns specify the orbits 
on which the 1- and 2-cycles evolve. The second column in- 
dicates the ty pe of 2-cycle using the classification scheme of 
section [I D 2j. The dynamical cycle that produces a signal at 



L ~ 6.09 has a prefactor of 2 with its 2-cycle class indicator to 
denote that it is the first repetition of a type-1 2-cycle. (In this 
case, each particle involved in the 2-cycle traverses one and 
one-half of the primitive orbit 7° before particle exchange.) 



orbit 7 = i(*2a). The particle of the 1-cycle completes 
one full motion on the orbit while the particles of the 2- 
cycle each traverse half the orbit and are then exchanged. 
The second peak (L w 4.17) occurs because of a dynami- 
cal cycle where the 1- and 2-cycles are on primitve orbits 
7° = ;j(*2a) and 7^ = ^(*4b), respectively. For the third 
peak (L 4.50), the 1-cycle is as in the first case, except 
the 2-cycle is type-0. Thus, as the particle of the 1-cycle 
completes one full motion on 7 = ^(*2a), the particles 
of the 2-cycle each traverse the full orbit and are then ex- 
changed. This is summarised in Table III where some of 
the dynamical cycles in this class are listed. In each case, 
the energies are divided according to Eq. (^). For L < 7, 
there are a total of 27 dynamical cycles. All structures 
in L-space can be accounted for in a similar manner and 
can be checked systematically by noting that due to the 
energy division between the particles, we expect peaks 

n%LP ) 2 /2. The 2-cycles are 



at positions L = y 

type-0 and type-1 for even (n+ ) and odd (n~ 2 ) integer 
repetition indices, respectively. 

The discrepancies between quantum and semiclassi- 
cal results are due to the problematic orbits mentioned 
above. The structure at L ss 5 is poorly reproduced due 
to the hetero-cycles involving 1-cycles that are station- 
ary and 2-cycles that are dynamical and evolving on the 
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5.6551 


i(*2a) 

2 V ^"7 


4a 


5.6560 


§(*2a) 


iriob) 


6.6031 


if 4b) 


±(*8b) 


6.8237 


i(*4b) 


4a 


6.8245 


|(*4b) 


!(*10b) 


6.9217 


i(*8b) 


|f2a) 



TABLE IV: Some dynamical cycles of the class (1, 1,0) that 
are responsible for numerical discrepancies. The first column 
gives the position of the signal in L-space arising from the 
dynamical cycle while the second and third columns specify 
the primitive orbits on which the 1- and 2-cycles evolve. 
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FIG. 10: Length spectrum for the class (0,0,1). Quantum 
(solid line) and semiclassical (dashed-dotted line) results for 
L < 6.25. Each peak is due to a dynamical 3-cycle of the full 
phase space. 



problematic orbits 7 = 4a and 7 = i(*10b). (The 
length spectrum of the hetero-cycles is shown in Fig. 
All other discrepancies are due to purely dynamical cy- 
cles and these are summarised in Table IV. 
We next consider the class (0,0,1). 



We numerically 
compare the Fourier transform of the dynamical 3-cycles 

F(S M (L)=HpU,i)(k)} (80) 
with its quantum analogue 

*W>,i)( L ) = ^{P(o,o,i)(*0 - P ( o,o,i)(fc)}- (81) 



1.50 
2.67 
3.00 
3.42 
3.80 
3.85 
4.50 
5.33 
5.52 
5.99 
6.00 
6.05 



Class 



1 
1 
2 
1 
1 
1 

2 
1 
1 

2(1) 
1 



i(*2a) 
i(*4b) 
§(*2a) 
i(*6b) 

3 a 
i(*8b) 
i(*2a) 
§(*4b) 
if 8c) 

5 a 
i(*2a) 
±(*10h) 



TABLE V: Some dynamical 3-cycles of the three-particle car- 
dioid billiard. The first column gives the position of the 
peak in L-space arising from the dynamical 3-cycle while the 
third column specifies the primitive orbit on which the 3-cycle 
evolves. The second column indicates the type of 3-cycle us- 
ing the classification scheme of Fig. ^. The dynamical 3-cycle 
that produces a signal at L = 6.00 has a prefactor of 2 with 
its class indicator to denote that it is the first repetition of a 
type-1 3-cycle. This situation is described further in the text. 



We use the first 1000 energies of /O( 0j o,i)(fc)- The smooth 
term is computed from Eq. ( |67| ) and trace formulas are 
computed using geometrical orbits with length L < 11. 
The result is shown in Fig. [l^. We now identify some 
of the peak structures with one or several of the orbits 
listed in Fig. |[ 

The first peak (L w 1.5) can be identified with a type-1 
dynamical 3-cycle consisting of all three particles evolv- 
ing on the orbit 7 = ^(*2a) with the same energy and 
exactly T 7 /3 out of phase and each completing one-third 
of the full motion on the orbit and finally being permuted 
as specified by r = (acb). The third peak (L sa 3) is due 
to a type- 2 dynamical 3-cycle where all three particles 
evolve on the orbit 7 = ^(*2a) as above, except that 
each particle completes two-thirds of the full motion on 
the orbit before being exchanged according tor= (abc). 
The peak at L « 4.5 is from a type-0 dynamical 3-cycle 
consisting of all three particles starting and evolving to- 
gether in phase on 7 = i(*2a), but each completing 
one full motion on the orbit and then being trivially ex- 
changed as presribed by any group element r G (0,0, 1). 
As an example of a higher multiple cycle, consider the 
first repetition of the type-1 cycle mentioned above. It 
is the same as before except that each particle completes 
one and one-third of the motion on the orbit before being 
permuted. This is summarised in Table |v| where some of 
the dynamical 3-cycles are listed. (The peak at L « 4.25 
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is completely undetected by the trace formula since it 
arises from a diffractive orbit. Such orbits require a sep- 
arate analysis since they are not included in the standard 
Gutzwiller theory As before, the discrepancy that 

occurs at L ss 5 arises from the two orbits 7 = 4a and 
7 ° = |(*10b). 

All structures in L-space can be accounted for in a 
similar manner. As a systematic check, recall that each 
dynamical 3-cycle can be mapped one-to-one with a pe- 
riodic orbit 7 of the one-particle phase space. If the orbit 
has length L 1 — n 7 L°, where n 7 is a repetition index, 
it is mapped to a 3-cycle where each particle executes a 
fraction n 7 /3 of the full motion on 7. We can then write 
n 7 /3 = i + j/3, where i = int(n 7 /3) (i.e. integer part of 
n 7 /3) and j/3 (j = 0,1,2) is the remainder. If j 7^ 0, 
then the orbit with length L 7 is associated with the ith 
repetition of a type-j dynamical 3-cycle. If j = 0, then 
it is the (i — l)th repetition. To determine peak posi- 
tions, we recall that all particles of the 3-cycle have the 
same energy E/3 and thus we expect peaks at lengths 
L = (Recall that a billiard orbit with length 

L 7 has action S 1 (£)/Ti = y/aeLj.) 



VI. CONCLUSION 

We began with the case of two noninteracting identi- 
cal particles evolving dynamically on periodic orbits and 
explained how the time-translational symmetry leads to 
families of periodic orbits in the full phase space. Us- 
ing the trace formula for continuous symmetries pc| , we 
obtained a trace formula for the two-particle resolvent 
consistent with the dynamical term of the semiclassical 
two-particle density of states jn). We also proved iden- 
tities for the s ymm etry-reduced densities of states (see 
also Appendix |A l|) w hich were stated without proof in 
our previous work |31] . Dynamical pseudo-periodic orbits 
(DPPOs) were defined and it was shown that their con- 
tribution to the semiclassical exchange term has the form 
of a one-particle trace formula. We also introduced two- 
particle heterogeneous periodic orbits in the full phase 
space. We discussed how the structure of these orbits is 
different in billiards and in analytic potentials and the 
explicit contribution of such orbits to the two-particle 
resolvent was determined. 

We have also demonstrated that the approach used in 
this paper yields results that are consistent with those of 
the convolution method |3l[]. In the convolution picture, 
one is faced with the asymptotic analysis of many convo- 
lution integrals and the further issue of spurious contribu- 
tions from them. The full phase space formalism, on the 
other hand, more easily generalises to arbitrary particle 
numbers. It is also more illuminating since it reveals the 
underlying structure of the periodic orbit families. (One 
useful property of the convolution approach occurs for 
billiards where significant higher-order corrections from 
hetero-orbits can be explicitly calculated (cf. Appendix 
). Most importantly, the convolution formalism does 



not accomodate interactions and it is necessary to use 
the full phase space if interactions are to be included. 
The symmetric/antisymmetric resolvent for the case of 
noninteracting particles can be expressed as a sum of re- 
solvents, one for each element of the permutation group. 
As shown above, trace formulas for the oscillatory compo- 
nents can be written as products over cycles, where each 
cycle is assigned to a periodic orbit of the one-particle 
phase space. The results of the paper were applied to a 
specific problem: three noninteracting identical particles 
in a two-dimensional cardioid billiard. We found that our 
semiclassical analysis correctly reproduced the quantum 
results and we explained how these results could be un- 
derstood in terms of classical structures in the full phase 
space. 

We have assumed that the single-particle dynamics is 
free of any continuous symmetry. If there are additional 
symmetries, then they must also be properly accounted 
for in the theory. If we restrict ourselves to the nonin- 
teracting case, then essentially the only difference from 
what we have presented above is that J and © become 
higher dimensional. (An example would be two particles 
in a disk billiard. In this case, there are four independent 
constants of the motion, any two of {E, E a , E^} and any 
two of {L z , L Za , L Zb }. Thus, the periodic orbits occur 
in 3-parameter families.) One future goal is to consider 
separately the important zeroth-order problem of har- 
monic oscillator potentials. The harmonic oscillator has 
a higher degree of symmetry than we are accounting for 
here (SU (d) in d dimensions). This project would require 
using the theory of Ref. |3(| which derives a trace formula 
for systems with more general symmetries including non- 
Abelian cases. 

As mentioned above, a major advantage of the for- 
malism presented here (as compared to the convolution 
formalism of Ref. |3lJ) is that it can be extended to in- 
clude interactions. Any interaction between the parti- 
cles destroys the periodic orbit families described above 
and replaces them with a discrete set of isolated orbits. 
For weak interactions, we can use a perturbative method 
p5| which is applicable to any situation where contin- 
uous symmetries are broken. This program will be ex- 
plored in a future publication |]32f . One could also study 
(zero-range) point interactions. Such interactions are of- 
ten considered as corrections to mean-field approxima- 
tions. Semiclassically, point-interactions can be under- 
stood using the formalism of diffractive orbits [^6|. One 
can imagine that such an interaction leaves the periodic 
orbit families (of the noninteracting system) largely un- 
changed, but introduces qualitatively new diffractive or- 
bits. We plan to explore this scenario in future work. 



Acknowledgments 

We thank Rajat Bhaduri, Stephen Creagh, Randy Du- 
mont and David Goodings for useful discussions. This 
work was financially supported by the Natural Sciences 



23 



and Engineering Research Council of Canada (NSERC). 



1. Symmetrised 2-particle Thomas-Fermi term 



APPENDIX A: TWO-PARTICLE 
THOMAS-FERMI CALCULATION 

We first discuss the smooth two-particle density of states 
and its decomposition into bosonic and fermionic densi- 
ties. Using the identity S(E — h(z a ) — h(zb)) = / de5(e — 
h(z a j)S(E — e — h(zb)), we can show that the leading- 
order smooth term for the two-particle density of states 
is the autoconvolution of the leading-order smooth term 
of the one-particle density of states (|l^) . We could verify 
this term- by-term in the expansion of p2{E), but we can 
do it more efficiently for all terms as follows. 

We use the partition function Z(j3) = Tr(exp (—(3H)) 
which is the Laplace transform of the density of states. It 
is convenient to work with the Wigner transform which 
is defined for an arbitrary operator A as 



A w (z) = J dx(q+ - 



A 



q -2 



exp 



in terms of which the trace is 



Tr{i} 



1 



{2irh) r 



dzAw(z). 



(A2) 



The trace of a product of two (but not more) operators 
is given by 



Tr{AB} 



1 



dzA w (z)B w (z). (A3) 



The Wigner transform of the evolution operator, 
exp (— (3H) w (z), can be written as an asymptotic expan- 
sion in powers of h, the first few terms of which are typ- 
ically retained and used as the smooth approximation 
to the partition function. Taking the inverse Laplace 
transform then gives the smooth density of states. In 
particular, the leading-order term of exp (— f3H) w (z) 
is exp (— f3H\v(z)), where the Wigner transform of the 
quantum Hamiltonian H^{z) is simply the classical 
Hamiltonian which we have denoted by H(z). (There 
are corrections to this if the Hamiltonian is not of the 
kinetic plus potential form.) The inverse Laplace trans- 
form of this expression yields the leading-order smooth 
term @. 

For two independent particles, the full quantum Hamil- 
tonian is the sum of one-particle Hamiltonians and since 
these are functions of independent phase space variables, 

exp (-PH) w (z) = exp (-/3/i) w (z )exp (-(3h) w (z b ). 

.(A4) 

Thus, the smoothed two-particle partition function is 
simply the product of smooth one-particle partition func- 
tions. By the Laplace convolution theorem, this implies 
that the smoothed two-particle density of states is the 
autoconvolution of the smoothed one-particle density of 
states. This same argument can be made for the exact 
density of states as an alternate proof of Eq. (||) . 



The bosonic and fermionic partition functions are 
Z±{P) = Tr(P± exp (— (3H)), where P± are the projec- 
tion operators defined in Eq. (]?]). The leading-order term 
is just the two-particle partition funtion. The next term 
Tr(U exp {—(3H)) requires slightl y m ore analysis and can 
be evaluated directly using Eq. (A3). We begin by find- 
ing Uw(z). 

It is shown in Ref. jl0| that for a one-particle sys- 
tem with a symmetry axis through the coordinate g, the 
Wigner transform of the reflection operator is 



(*) 



K) (z)=«h5(q)8(j>), 
w 



(A5) 



where p is the momentum conjugate to q. We map our 
problem onto that one as follows. First, suppose that 
the one-particle system is one-dimensional and define the 
Jacobi coordinates: q = q a — p = (p a — pb)/2, Q = 
(q a + q b )/2 and P = p a + Pb- Then, exchanging a and 
b is equivalent to reflecting in q so that the variable q 
in the above equation is replaced by q a — qb and p is 
replaced by the conjugate momentum (p a —pb)/2. Then, 
U~w(z) = 2irhS(q a — qb)S(p a — pb). If the one-particle 
system is higher dimensional, then U is the product of 
one such inversion in every component. All of them are 
independent so that the final result is the product of the 
individual ones (for the same reason that Eq. (A4) is 
multiplicative). The final result is 



U w (z) = (2Trh) d S{z a ~Zb), 



(A6) 



where the delta function represents the product of all 2d 
delta functions (two for each component). Equivalent 
results can be found in Ref. jl7| . 
Then, 



Tr (Uexp(-/3H) 
1 



(2TTh) 2d 
1 

{2TTh) d 



dzC/ w (z)exp (-/3H) w (z) 
dz a exp (-2f)h) w (z a ), 



(A7) 



where we 
Eq. 



used the delta functions from Uyy{z) in 



(A6) to do the integrals over the Zb variables and 
the mult iplicative property of the W ign er functions as in 
Eq. (A4). The first few terms of Eq. ( A7) give the smooth 
approximation to the one-particle partition function eval- 
uated at 2/3. Under the inverse Laplace transform, this 
becomes p\{E/2)/2 and we conclude 



P ± {E) = \Ue)±\p, (f 



(A8) 



consistent with (|Io|). 
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APPENDIX B: STABILITY MATRIX OF 
PSEUDO-PERIODIC ORBITS 

We first prove that det(M y - 7) = 4det(M 7 - I), 
where My is the stability matrix of a pseudo-periodic 
orbit 7' of the full two-particle phase space and M 7 is 
the stability matrix of the corresponding periodic orbit 
7 in the one-particle phase space. This admits various 
generalisations which are used in the main discussion. 
The (type-1) dynamical pseudo-periodic orbit (DPPO) 
consists of both particles evolving for half of the single- 
particle period T 7 /2 followed by the symplectic mapping 
u that exchanges the two particles. 

We define coordinates as follows (cf . Fig. |ll|) . For par- 
ticle a, we define an initial section S a such that the phase 
space flow is transverse to it and all points on the sec- 
tion are at equal energy. We define a coordinate pointing 
along the orbit which we call rj a . Without loss of general- 
ity, we can take dr\ a jdt = 1. We also define a coordinate 
transverse to the constant h a surface (but in the phase 
space of particle a) which we call re a . If we consider it 
to take the values of h a , then it is canonically conjugate 
to r\ a and has zero time derivative under the flow since 
h a is conserved. The remaining (2d — 2) coordinates for 
particle a lie on the section S a and will be collectively 
denoted by £ . As the flow evolves, changes in the £ a 
coordinates are described by the (2d — 2) x (2d — 2) sym- 
plectic stability matrix (for the one-particle dynamics) 
N a . Similarly, we define r)b , Kb, £b and iVf, for parti- 
cle b. We also need a way to connect coordinates on S a 
to those on Ej,; we will take them to be such that they 
are connected by parallel transport so that, for example, 
the stable and unstable manifolds are mapped onto each 
other. 

We start by defining the symplectic transformation 



Va + Vb 



77 = 



v = r)a- Vb, 



K — K a -f- Kb, 
j. K a Kb 



(Bl) 



The monodromy matrix My describes the linearised mo- 
tion of small perturbations around a DPPO 7' of the 
full phase space. In particular, if Y = (k, 77, v, f, £ a , 
then <5Y = My 5Yo. Consider an initial slight change 
in 77 by the amount St]q while keeping all other coordi- 
nates constant. This implies that both r\ a and r/b in- 
crease by Stjq. After time evolution for T 7 /2 and par- 
ticle exchange, 8rj = 6r]o while all other coordinates are 
unchanged (in particular, the transverse coordinates are 
unaffected) . Now consider an initial small change in k by 
the amount Skq. This implies that K a and Kb change by 
5kq/2. After integrating for time T 7 /2 and interchanging 
the particles, we observe that 5k = Skq. However, this 
change of value in k does affect the 77 coordinate. Under 
this change, the period of the orbit 7 also changes; let 
T 7 denote the derivative of this period with respect to 
the single-particle energy. Since we are only integrating 
for half of the period, and the single-particle energies are 




FIG. 11: The coordinates of particles a and d on a (type- 
1) dynamical pseudo-periodic orbit 7' of the full phase space 
(which can be mapped one-to-one to an orbit 7 of the one- 
particle phase space). E a denotes an initial section for particle 
a, rja is the coordinate along the orbit (the coordinate trans- 
verse to the h a surface denoted by K a is not shown) and £ a 
are the (2d — 2) remaining coordinates for particle a which lie 
on E a . (All points on E a are at equal energy.) Similarly, for 
particle b. 



changed by Skq/2, we find that Sr/ = —T^&kq/A. (The 
minus sign indicates that if the period increases and we 
integrate for the same amount of time as before, then 
the particles will fail to execute a complete loop, corre- 
sponding to a negative value of 77.) Thus, the monodromy 
matrix of the pseudo-periodic orbit 7' in the full phase 
space has the form 



My 



( 



\ 



1 









1 









4 









M 7 



(B2) 



We are interested in calculating det(M 7 ' — I) to eval- 
uate the Gutzwiller amplitude. Note that the matrix 
Afy involves only the (Ad — 2) phase space coordinates 
other than 77 and re. We can understand the calculation 
up to now as follows. The transformation (Bl) can be 



thought of as a transformation to center-of-mass coordi- 
nates. We have removed the center-of-mass coordinates 
77 and re from consideration and are only left with the 
relative coordinates v and £ (as well as all the transverse 
coordinates £ a and £&.) It is reasonable that only the 
relative coordinates are important for determining the 
stability. 

The next two coordinates we consider are v and £• 
Let us start with v. A small initial change in v by the 
amount 5vq implies that 77 Q changes by Svq/2 while 77b 
changes by —Svq/2. After integrating for time T 7 /2, this 
remains unchanged, but after particle exchange, the fi- 
nal values of <577 a and 5rjb are changed in sign so that 
the corresponding diagonal matrix element of My is — 1 . 
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Similarly, the diagonal matrix element corresponding to 
the C coordinate is also —1. As before, an infinitesmal 
change in £ implies an infinitesmal change in u. In this 



case, the corresponding matrix element is T^. 
we can write 



M. 



7' 



V 



-1 



V 

7 



Therefore, 



(B3) 







N 



Then, det (My — 7) = 4det(/V — I), where we use ap- 
propriately dimensioned identity matrices on each side 
of the equality. It remains to calculate the determinant 
of the (Ad - 4) x (Ad - 4) matrix N — I. The matrix TV 
involves only the coordinates £ a and £&. Since these two 
sets of coordinates live on different sections, we cannot 
immediately define a mapping between them. To do so, 
we note that we have defined coordinates on the two sec- 
tions so that the exchange operation is a simple mapping 
of the form 



E 




(B4) 



where I is a (2c?— 2) x (2d— 2) identity matrix. In terms of 
these coordinates, the one-particle stability matrices N a 
and Nb are such that N a N b — M 7 , which is the stability 
matrix of the full periodic orbit for the one-particle dy- 
namics. The combined operations of flow and exchange 
give 



and 



N = 



det(N - I) 



N b 

N„. 



(B5) 



det 



det 



N a -I 

N a -I 
-I N b 



det (M 7 - 1), 



(B6) 



where after the second equality, we interchanged rows to 
put the matrix into a more useful form. (The matrix 
has even dimension so there is no sign introduced as a 
result of this interchange.) The final equality of Eq. (B6) 



requires the following identity. If a matrix C has the form 



C 



then det(C) = det (AS - I). 
tiplying C by the matrix 



C' = 




(B7) 



This can be shown by mul- 




(B8) 



After multiplying them together, the product is block- 
diagonal with (AB — I) in one block and (BA — I) in 
the other. These have equal determinants. Since C has 
the same determinant as C, it follows that [det(C)] 2 = 
[det(AB-J)] 2 and thus we have identified the two deter- 
minants within a sign. The sign follows from the fact that 
the contribution to the determinant of the fully diagonal 
term ^[^4^73^ should be positive. Thus, we conclude 
that det(My - I) = 4det(M 7 - I). 

It is a straightforward extension to generalise this re- 
sult to a cycle with n particles on an orbit. We first 
have to find some appropriate set of variables so that we 
may isolate an upper-lef t bl ock of the monodromy ma- 
trix in analogy to Eq. (B3). This comes from the 2n 
coordinates i] and k. As in the previous case, we sepa- 
rate the variables into center-of-mass coordinates (which 
subsequently play no role) and a set of relative (Jacobi) 
coordinates. Using similar arguments to the n = 2 
case, the determinant of the upper-left block is then 
n 2 . The contribution from the rest of the matrix (i.e. 
the lower-right block) comes from the transverse coor- 
dinates. In terms of these transverse coordinates, the 
single-particle stability matrices N a , N b , . . . , N n are such 
that N a N b ---N n = Af 7 , which is the stability matrix 
of the full periodic orbit for the one-particle dynamics. 
Through a sequence of manipulations and transpositions 
similar to the n = 2 case, the determinant of this lower- 
right block (i.e. det(7V — /)) can then be reduced to the 
form 



det 



f Na -I 
N b -I 


\ -I 



\ 


Nn-i -I 
N n J 



(B9) 



which is a generalisation of the n = 2 case. It can 
be shown that det (TV — I) = det(M 7 — I) and thus, 
det (My — I) — n 2 det(M 7 — I). If the orbit is not prim- 
itive, but is a repetition of some simpler orbit, then we 
can absorb this into the definitions of the single-particle 
stability matrices and carry through all of the manipula- 
tions as before. The result is unchanged. 



APPENDIX C: MONODROMY MATRIX AND 
PHASE INDICES OF A HARMONIC 
OSCILLATOR 

It is shown in Ref. § that the monodromy matrix for 
a primitive orbit along the x-axis is 



M 



cos(u y T) ± sm(uj y T) 
-ojy sin(w a T) cos(u> y T) 



(CI) 



where u> y is the frequency of the y-motion and T — 2tt/uj x 
is the period of the x-motion. This is derived by integrat- 
ing the harmonic oscillator equations of motion for time 
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T. Then, det(M - I) = 2 sin(w y T/2). In higher dimen- 
sions, the monodromy matrix is simply block diagonal so 
that det(M - /) = flj 2sin(wjT/2), where the product 
is over the directions other than x. 



The only role played by the x variable above was to 
specify the time of evolution in the determination of the 
arguments of the sinusoids. It was not important that it 
be a harmonic motion, it is enough that it be periodic. 
These results apply for any periodic orbit with period 
T as long as the transverse motion is harmonic. This 
exactly describes a heterogeneous orbit. This justifies 
the amplitude factor in the denominator of Eq. (^8|) for 
d = 1. In higher dimensions, the stability is given by 
both the single-particle monodromy matrix of the evolv- 
ing particle and by the harmonic motion of the stationary 
particle about its potential minimum. But these motions 
are uncoupled and so are simply multiplicative in their 
combined contribution to the amplitude. 



The phase factor can be determined in an analogous 
manner from the exact harmonic oscillator trace formula. 
For the primitive orbit along the x-axis, the phase index 
is 3. A factor of 2 arises from the two turning points expe- 
rienced by the periodic orbit in traversing the ^-motion 
independent of the harmonic motion transverse to the 
orbit. The remaining factor of 1 can be attributed to 
the harmonic y-motion and is related to the sign of the 
determinant of the monodromy matrix. For heteroge- 
neous orbits, this means that we should simply include a 
phase factor of — 7r/2 for the transverse harmonic motion 
in addition to any phase factors from the single-particle 
motion along the periodic orbit. In higher dimensions, 
each transverse direction is independent and the phase 
index is additive. This accounts for the phase factor of 
—dir/2 in Eq. (28). The fact that each transverse direc- 
tion is uncoupled from all the rest as well as from the 
single-particle dynamics transverse to the periodic orbit 
allows us to simply multiply the amplitudes and add the 
phase factors. 



Finally, if the potential is a local maximum in one of 
the directions, this corresponds to the case of an unsta- 
ble harmonic oscillator. It is straightforward to show that 
its contribution to det(M — I) is 2 sinh(wj / T/2). Further- 
more, its phase index is trivially zero since an unstable 
periodic orbit running along a ridge does not fold back 
on itself and introduces no caustics in phase space. This 
fact is also consistent with the trace formula for an unsta- 
ble harmonic oscillator as described in Ref . |6| . In higher 
dimensions with a mixture of stable and unstable direc- 
tions, we continue to multiply the amplitude factors and 
add the phase indices of the separate directions. This 
fully accounts for the modifications described immedi- 
ately below Eq. @. 



APPENDIX D: SYMMETRISED iV-PARTICLE 
THOMAS-FERMI TERM 

We now discuss the smooth contribution to Tr(£/ r G) 
and we sh all e valuate it using Wigner transforms as in 
Appendix A 1 . We need to determine the smooth ap- 



proximation to the symmetric/antisymmetric partition 
function 



Z±(/3) = Tr(^P ± e X p(- ) SF) 

= 4rE( ±:L ) nTTr (^ ex P(-^) 



Nl 

T 

(Dl) 

Since each group element can be decomposed into inde- 
pendent cycles, 



(£t)w(z)=I]>( 



Zk 



(D2) 



where k indicates the different cycles comprising the 
group element r, fk is a function we discuss below and 
Zk denotes the phase space coordinates of the n& parti- 
cles being permuted by that cycle. (For each group ele- 
ment, the unique decomposition into cycles also provides 
a unique decomposition of the phase space into the sub- 
spaces corresponding to the cycles.) The function fk can 
be specified without loss of generality by choosing to label 
the particles being permuted by the cycle as 1, 2, . . . , rik 
(i.e. 1 — ► 2, 2 — > 3, . . . , rife — > 1) and to leading order in 

fk(zk) « (2Trh) ( - n x-V d 5(z 1 ~z 2 )S(z 2 -z 3 )---6(z nk _ 1 -z nk ). 

_ (D3) 

The first group element of the sum in Eq. (Dlh is the 
identity element for which the decomposition into cycles 
is the trivial one where each particle is in a cycle by itself 
so that all of the fk are identically unity. Integrating the 
smooth approximation to the Wigner function of e~@ H 
yields the smooth iV-particle partition function. Us- 
ing the generalisation of property ( |A4| ) , we observe that 
the leading-order term of Z±(j3) is just the Nth. power 
of the single-particle smooth partition function Zi((3) N 
and under the inverse Laplace transform, this is just the 
(N — l)-fold convolution integral of the single-particle 
smooth density of states. The prefactor of 1/Nl comes 
from the projection operator (JO) and we conclude that 
the identity term is 0(1/N\K ). The first correction 
will come from group elements that consist of one 2-cycle 
and (N — 2) 1-cycles. The contribution from this class 
will have the form Zi((3) N - 2 Zi{2[3). Compared to the 
leading-order term, this class contributes to the density 
of states with relative order 0(Nh d /2). The factor of N 
is due to the fact that this class has N(N—l)/2 members 
and they all contribute identically. The factor of 2 comes 
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from the inverse Laplace transform since the argument of 
one of the single-particle partition functions is 2(3. The 
general structure then emerges. For an arbitrary group 
element, the contribution to the smooth partition func- 
tion is Y\ k Zi(rikf3). It contributes to the smooth den- 
sity of states with relative order 0{Tv- N ~ mT ' w T ), where 
m T is the number of independent cycles in the decom- 
position of r. The factor w T is the size of the class (a 
combinatoric factor which can be found from Eq. (1-27) 
of Ref . |34| ) divided by a factor arising from the inverse 
Laplace transform which equals Ylk n k- 

As a formal expansion in powers of h, this may be in- 
consistent since some of the neglected corrections from 
the first few group elements may be of more significant 
order than the leading-order contributions of later group 
elements. However, for large N, one could easily imagine 
that the combinatoric factor w T offsets this effect. Keep- 
ing the leading-order h term of all the group elements 
then guarantees that one has a good approximation re- 
gardless of the relative size of 1/h and N . 



APPENDIX E: HIGHER-ORDER h 
CONTRIBUTIONS FROM HETERO-ORBITS IN 
2-D BILLIARDS 

In this paper, we only discuss leading-order contri- 
butions to the oscillatory part of the density of states. 
For billiards, hetero-orbit families generally have higher 
dimensionality than dynamical orbit families and the 
corrections from the former can be quite significant. 
The formalism of Ref. |29] cannot be used to determine 
these corrections. We calculate correction terms arising 
from hetero-orbits using the convolution method as in 
Ref. |HJ. The first few terms of the asymptotic scries 
can be determined by convolving the Weyl expansion Jl3| ) 
term-by-term with a two-particle trace formula. As a for- 
mal expansion in powers of h, this is inconsistent since we 
do not include corrections to the one-particle trace for- 
mula ( |l5| ) itself. However, our numerics seem to indicate 
that the corrections to the Gutzwiller trace formula are 
negligible in this case, otherwise we could not reproduce 
the quantum results with such accuracy. 

The contribution from the first type of hetero-orbit 
where one particle evolves while the others are stationary 
is calculated from 

pf(E) = p 1 {E)*{p 1 *p 1 ){E) 

pi(e)(pi*pi)(.E-e)d£. (El) 

After convolving Eq. (|l^) with the oscillatory function 
(pi * Pi){E) (which has been calculated in Ref. fsifl), we 
find there are nine integrals to do, but three of these 
are trivial because of a delta function in the integrand. 
The remaining six integrals require careful analysis. As 
an example, we obtain the asymptotic expansion of one 
integral. The others are calculated in the same manner, 



but we forego the details. The first correction to the 
leading-order result ( f49[) comes from two terms; the first 
(second) is the convolution of the area (perimeter) term 
of p~i(E) with the perimeter (area) term of (pi * p~i)(E). 
The integral involved in the first term is 



Iac(E)=J (E-s)- 1 ^ cos (aVE~e + b^)de 7 (E2) 

where a = y / aL 7 , b = — ct 7 7t/2— 7r/4. We want to perform 
a local analysis about e = 0. The reason for this is that 
small e corresponds physically to the situation where the 
stationary particles have little energy and most of the 
energy belongs to the dynamically evolving particle. For 
e w E, we have the opposite situation which does not 
make sense physically. We first make several changes of 
variable to simplify the calculation. A first change of 
variable u = (E — e) 1 / 2 removes the square root in the 
argument of the sinusoid. Next, we wish to make a Taylor 
expansion about the point u = yE and to facilitate this, 
we make a second change of variable x = \fE — u. The 
integral then becomes 

rVE 

Iac{E) = 2 J (VE - x) 1/2 cos (a(VE~ - x) + bj dx 
«2Re|e wv/¥+b ^ {s/E - xf^e'^dx^ . 



(E3) 



The integral J E (■ ■ ■ ) = / °°(- • ■ ) — //gO • • ); the second 
integral is an endpoint correction, but asymptotic in E, 
this correction is negligible. Thus, we are justified in 
replacing y/~E with oo in the second line above. At this 
point, we obtain the Taylor expansion of f(x) = {^/E — 
x) 1 / 2 about x = 0: 



fix) = E 



-E-VS 



.£-3/4 



X* + . . . . 



(E4) 



Since the final correction [Ikc(E)] in the expansion (El) 
is 0(E^ 1 ^ 4 ), it is only necessary to include terms in the 
Taylor series to 0{E Ll / A ). Thus, 



{y/E-x) l ' 2 e 



*c\x 



lo 

i 



c dx (E5) 



E 



1/4 



a 



-E 



-1/4 



and asymptotically, 
2 



Iac(E) ~ - I £ 1/4 cos (aVE + b - J 
a L V 2 

£-1/4 



2a 



cos ( aVE + b 



(E6) 
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An equivalent approach is to evaluate the integral ex- 
actly and then replace the resulting functions with their 
asymptotic forms. Evaluating the integral (E2) at the 
upper limit using this method then corresponds to the 
situation where one of the stationary particles has all of 
the energy while the dynamically evolving particle has no 
energy. Physically, this does not make sense and this is 
evident mathematically since the contribution that comes 
from evaluating the integral at the upper endpoint is a 
smooth function of E and is therefore spurious in the 
sense discussed in Refs. [pll |37t|. The result is only mean- 



ingful if we drop this smooth contribution. This is justi- 
fied since we know that any smooth contribution to the 
density of states is already contained in the ps{E) term. 
This is completely equivalent to what is done above (the 
spurious smooth contribution above is the endpoin t co r- 
rection that was dropped in the second line of Eq. (E3)). 
All six convolution integrals can be analysed in this way. 
Collecting the contributions from all six integrals, the 
expansion up to 0(1/H 3 ^ 2 ) is 











det (M 7 


-') 


( ^ 


C 2 




f 32 



a 3/2^2ffl/2 
8^ 3 L2 



cos ($) 



■AK 



$ _ 



a^ACE 1 / 4 

p — ; — r/9 cos 

^CE- 1 / 4 ( 3A 



3tt 

T 



(E7) 



/C cos 



(•-3 



where $ = \fotEL 1 — er 7 7r/2. The contribution from the 
second type of hetero-orbit where one particle is station- 
ary while the others evolve is calculated from 

pf{E) = ME)*(pi*Pi)(E) 

f ^ 

= / p 1 {e)(p 1 *p 1 ){E-e)de. (E8) 
Jo 

I 



After convolving Eq. ( |l3| ) with the formula for {p\*p~\) (E) 
(which has been calculated in Ref. |3l| using stationary 
phase asymptotics), we find there are only two convolu- 
tion integrals which require analysis. These are evaluated 
asymptotically using the same technique as above. The 
final result (up to 0(l/h 3/2 )) is 



pf{E) = £ 



rO rO ( T 2 , r2 



,-1/4 



Ti.72 ^|det(Af 7l -I) \l det(M 72 -I) 
aC 



a^AE 1 / 4 / 3tt 



(E9) 



n\ a^E- 1 ^ ( A 
2 ) + - 



r 



where L\2 



L\ + L 72 , CT12 = o- 71 + cr 72 and $ i2 = VaEL 12 - a 12 7r/2. 
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The division by TV of the energy argument simply states 
that the total energy must be evenly divided among all of 
the particles. The set of orbits corresponding to this cycle 
is clearly the same as the set of orbits of the one-particle 
dynamics (almost by definition) and the amplitudes and 
actions are the same as the one-particle case since the TV 
particles collectively execute one complete motion (or a 
multiple repetition) of the periodic orbit. This is observed 
for TV = 3 in section [V B . 

the operator P% that projects onto 
[ z/\G\) J2 g Xx(g)U$, where the sum 
is over all group elements g £ G, dx is the dimension 
of the irrep, \G\ is the order of the group, Xxid) i s the 
character of the group element g in the irrep T and U g 
is the operator that transforms ^ as prescribed by the 
group element g G G. Permutation operators are unitary 
(Sjv is a unitary grou p). For A± irreps, this reduces to the 
operator of Eq. (kll) and for the irrep £ of 53, using the 
information provided in Table §, P £ = l(2I-U T1 ~U T2 ), 
where n/a G (0, 0, 1) and I denotes the identity operator. 



For a discete group G 
the irrep 1 is Pi 



